{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b21d7059-0777-4afe-827a-624b5a950d27",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from linearmodels.panel import PanelOLS\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "26e19683-cf21-4109-8290-6c5aa065c347",
   "metadata": {},
   "outputs": [],
   "source": [
    "panel_df = pd.read_excel(r\"\\df_analysis_1.xlsx\")   # 복지패널"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0548ffe5-52d8-4fd2-b9e9-3cd778fe69ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "cross_df = pd.read_excel(r\"\\df_analysis_2.xlsx\")      # KGSS 횡단면"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a65da0a1-bd52-473a-8f21-ce19221437c4",
   "metadata": {},
   "source": [
    "# 1. Analysis 1 and 2 plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4cf57363-fcad-43f1-ab88-48a87757aa49",
   "metadata": {},
   "outputs": [],
   "source": [
    "def safe_log(s):\n",
    "    s = pd.to_numeric(s, errors=\"coerce\")\n",
    "    out = np.where(s > 0, np.log(s), np.nan)\n",
    "    return pd.Series(out, index=s.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "346b301c-96db-4ae2-8d9c-65ad5645fe93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---------- 2) panel prep (two-way FE) ----------\n",
    "# cast numerics\n",
    "panel_numeric_cols = [\"id\",\"year\",\"ideology\",\"divorce\",\"housing_price\",\"income_1\",\n",
    "                      \"homeownership\",\"education\",\"sex\",\"age\",\"religion\"]\n",
    "for c in panel_numeric_cols:\n",
    "    if c in panel_df.columns:\n",
    "        panel_df[c] = pd.to_numeric(panel_df[c], errors=\"coerce\")\n",
    "\n",
    "# safe logs\n",
    "panel_df[\"ln_house_price\"] = safe_log(panel_df[\"housing_price\"])\n",
    "panel_df[\"ln_income\"]      = safe_log(panel_df[\"income_1\"])\n",
    "\n",
    "# key vars must be present\n",
    "panel_df = panel_df.dropna(subset=[\"id\",\"year\",\"ideology\",\"divorce\"])\n",
    "panel_df = panel_df.set_index([\"id\",\"year\"]).sort_index()\n",
    "\n",
    "# design matrix: DO NOT add constant in FE\n",
    "X1_cols = [\"divorce\",\"ln_income\",\"homeownership\",\"ln_house_price\",\n",
    "           \"education\",\"sex\",\"age\",\"religion\"]\n",
    "X1 = panel_df[X1_cols]\n",
    "y1 = panel_df[\"ideology\"]\n",
    "\n",
    "# fit two-way FE, drop absorbed regressors\n",
    "res1 = PanelOLS(y1, X1, entity_effects=True, time_effects=True,\n",
    "                drop_absorbed=True).fit(cov_type=\"robust\")\n",
    "\n",
    "# divorce coefficient from FE\n",
    "if \"divorce\" in res1.params.index:\n",
    "    b1 = res1.params[\"divorce\"]\n",
    "    se1 = res1.std_errors[\"divorce\"]\n",
    "else:\n",
    "    # fallback: time FE 제거(개체 FE만)로 재추정하여 divorce 회수 시도\n",
    "    res1_alt = PanelOLS(y1, X1, entity_effects=True, time_effects=False,\n",
    "                        drop_absorbed=True).fit(cov_type=\"robust\")\n",
    "    b1 = res1_alt.params.get(\"divorce\", np.nan)\n",
    "    se1 = res1_alt.std_errors.get(\"divorce\", np.nan)\n",
    "\n",
    "# ---------- 3) cross-sectional prep (Logit) ----------\n",
    "# force numeric\n",
    "cross_numeric_cols = [\"conservative\",\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"]\n",
    "for c in cross_numeric_cols:\n",
    "    if c in cross_df.columns:\n",
    "        cross_df[c] = pd.to_numeric(cross_df[c], errors=\"coerce\")\n",
    "\n",
    "cross_df = cross_df.dropna(subset=[\"conservative\",\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"])\n",
    "\n",
    "y2 = cross_df[\"conservative\"].astype(int)\n",
    "\n",
    "# year dummies; ensure float dtype for statsmodels\n",
    "X2 = pd.get_dummies(cross_df[[\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"]],\n",
    "                    columns=[\"year\"], drop_first=True)\n",
    "X2 = X2.apply(pd.to_numeric, errors=\"coerce\").astype(float)\n",
    "X2 = sm.add_constant(X2, has_constant=\"add\")\n",
    "\n",
    "# Logit (show log-odds coefficients directly)\n",
    "logit_res = sm.Logit(y2.values.astype(float), X2.values.astype(float)).fit(disp=False)\n",
    "\n",
    "# divorce coefficient from Logit (log-odds)\n",
    "# find column position of 'divorce' in X2\n",
    "div_idx = list(X2.columns).index(\"divorce\")\n",
    "b2  = logit_res.params[div_idx]\n",
    "se2 = logit_res.bse[div_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3b46af7d-f9af-43bd-b24e-5001e83730bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABL8AAAGGCAYAAACaMSNVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAmEJJREFUeJzt3QWYVGX7x/Gb7hAJaWlEEVDBQMXAwlcRJcRAxe56LQzsxm7FDkxMbBS7QARBEJUSFAnphvO/fs++Z/6zs2dmZ3ZndmbPfj/XtQ5OnDlz8j738zz3Ked5nmcAAAAAAABACJXP9gwAAAAAAAAAmULyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFF8gsAAAAAAAChRfIrwI8//mgVKlSws88+Oy3T6927t9WuXdv+/vvvIn1+6623tnLlyuX7q1KlirVo0cIGDRpkn3/+uZU2s2bNcr9Dvy0Ve+21V4FlEft33nnn5ftMYe/3/z799NNibytPPvlk0t8X/afPIT02bdpkHTt2tJYtW9qaNWvSPv2rr766wPqrWrWqNWzY0Lp06WLHH3+8Pffcc7Z27doCn73qqqvc+3fZZZekvmv06NHu/fXq1QucXrboN8Yug4oVK1qDBg1sv/32s6effto8zyv29/j7k74vHceP4n42Xfxlliv842pRjoFh4W/THIuTj4ni7Z+F+eKLL9znLr744mLP2+LFi+2mm25y2/BWW21llStXdvHWdtttZyeffLKNHTu22N+BzPJjXB2bSxsdMzXv2v6K8ptz8XiTrXMksVV+v/zyi11wwQXWrVs323LLLa1SpUrucdddd7XLLrvMvY7cU9TzYqatWLHCatas6ebtvffeS+ozXbt2de+/9dZbLUwqZnsGcpECvGrVqtmVV16ZlundfPPN1r17d3eweuKJJ4o8nZ49e1rbtm3dv5cuXWo//PCDvfTSS/byyy/b7bff7g6SZYVOhNopg/To0SPw+QMOOMAFx/Ekei3ZbUXr57jjjgsM9n///Xdr06aN7b777gVe99crik8Xaddff70NGDDAHbCHDx+eke9p1KiRHXjggZGE27Jly2zatGn21FNPuT8lYe+991478sgjI5854YQT3Lx9++23NnXqVOvUqVPC73j88cfd49FHH+2CwFwTvT0rgPz555/to48+cn9vvPGGOz5pfZQkBeyzZ8+2mTNnZjXBhdTpAnLcuHH2ySefpHwxWVgwqn1Px+ZcvNgsazGRjhkHH3yw3X333S5B1a5duyJN55lnnrEzzjjDVq5c6RoEde5v2rSpa/TQsfixxx5zfzoX6FiEkqcLQJ0PFfvm2sUgcvP8WtZjq40bN9pFF11k99xzj23evNkl6HQNqcSXrv3Gjx9v33zzjYtvdQw966yzsj3LZYYSw61atXKN66UtWV+rVi13LnzyySfd9u/vY/FoO/vpp59cw/aQIUPCtbw85PPyyy+ru4J30UUXpXW6//nPf7xy5cp5EydOTPmzLVu2dPP0xBNP5Ht+zZo13pAhQ9xrFSpU8KZPn+6VFjNnznTzrd+Wil69ernPDR8+POnP6P36++STT7xsbSvHHXece68eUTI6d+7sVatWzfvrr7/SOl1te1qX2haD/Pbbb94xxxwT2e7uv//+fK/vu+++7vkLL7ww4fdovitWrOje++OPP3q5JNH2/MADD0R++8iRI4v1PTrmBX3P+vXrvV9++cUt63jHSx1jgiT6bEnxl0+umD17tlsmq1atyup8+Mf3dB+r421H0ebPn++WwdKlS9P63aVdovNcMss1nh9++MF9tl+/fkWarwcffNB9XnHVJZdc4i1btqzAe6ZMmeINGDDA69q1a5G+A+k7V8TGr9F0LNa+p2NzaaNjVaJ4INW4PhcU5/xaHMRWeQYNGuTmrXbt2m772LhxY77XN2/e7L3//vveDjvs4J177rlZm8+yKJlrV8UQ2n8UU+Sazz//3M1/lSpVvMWLFyd87xlnnOHe27dv3xK/1s80hj3GuPPOO93jiSeemNbpanq65lGWPl3UWnH//fdbjRo1XOvIa6+9lrZpI3vbCtJj6NChrgfAI488UuK9odQjQS13cu6559off/wRed3fXp599lnXwhePhg7qdXV5j9fLMRedfvrp1qtXL/fvTPW2UPd/DW3Vsi7Jz4aVhtBrmVSvXj3bs5I1jRs3dsugTp062Z6VMnGe23HHHV0PbvUQTbVFWL1AzjnnHPfvESNGuN71GuoYS70/dAxKZ9yF9NOxWPuejs3Ivlw9R5aF2Eo9cl588UW3Dj744APXWzK297yGoe2///6u95dK3yC3KIbQ/qOYIteo13WHDh1s3bp1bghxPHr9hRdeiFxLhQ3Jr5i6Fl999ZUbM66NI9aGDRvcQVXdZLVhK9jSUAC9V4HY/Pnz405bXfzr16/vNqYlS5akbZ41ftefVz+AXLhwoesu26dPH9fdUPOoed1pp53slltuiTu+PboOzauvvup2En1OyTUNuRwzZkzc+dCJRMMLNFRFXXQ1BEHfrQvhuXPnWlnbVpIxadIkt7y33377Aq+pO7O/PtSFO9qcOXPc87pgjfXnn3+6ISoaRqLkqA7CWncPP/ywS5CmaocddnDfNWHChHzP//PPP1a+fPm4dVv22Wcf91p0vZVUt0sNE9VJf4sttrDVq1fHncdtt93WfVfs9qn9VN119dsTBUKZcsMNN1iTJk3cd/sXkHL44Ye7fWTBggX2zjvvxP28P0S6NCZXdWErsRe1umjV8AR1gdYxQsth3333TTlJFlSTxK+zoCEZom0sqKZfYfVMtK3ddddd7vinbU/zqfk95JBD7Pnnn8/3Xn2Xtl1t79of9d66deu6z2q705CFdNK86UJf+6W6sCtZpe3/iiuusH///TfhctL+f8cdd7iA36/7kGzNr48//thttwrmVFNJNVj69etnX3/9dVrOJX7dHA15lL333jtuTUQNqdUxThctOqdqmTdr1sxdBHz//fcF5kW/XducaMhM9HSjh1YWVvNr1KhRblv1z2/aJhQU/vrrr4XWMdIwTl2saHvScU/rTxdgYTrPqfbWmWeeGdkPtHzOP//8wO0yeplrH3nwwQdTmi/tc4rHlDyLrfEZZM8994y7feo4qxo6OlfG1p367rvvbODAge447m/3Og58+OGHcS8YbrvtNnf80/6pz6icgoYs6TwZG/vNmDHDbUM6VmmZab/UclO8WJQSGdo3NH8aNqaLZ21vigWOOeYY++yzz9Kybyd7jPSPPdrnRPtg9L6n+k7J1PwqzjFPDc5q/NL60LFH61j7YbzfpvWt9aThs379OC1L/S4t22wrSnynZaCEiuIsLTsNnTvooIPcfh2vXllxzq8lIayxldaVfpvo2mnnnXdO+H7t4zp2xUr1uBV97tNw1mOPPdZt/9qvlXDUvqZjWywdu7V/aRtU3KP58euzaTsN2p9TvVaM3kZ1LFBtt2222cZty9o+FVPqdR2HEtVu0/av96mxxadrK5VF0fxruLyWlfYP1ekOikm1nDSvov0gtj5doppf2kf1XKKhhjqHanloPnS9FE3HOs2r4h7/ONi5c2c31DfR9VGQE/+33ftDf+PVxNN3ajvQNZtP57Bhw4a5Y7DmQfOi46uuWWNrLCe7vKKHWeq6zY8htH2oXFGi3EORZbvrWS656qqrXPe8K664IvD1uXPnutfr1Knj7bLLLq47fZ8+fbwmTZq45xs0aODNmDEj7vT79+/v3vf888+ntXt027Zt3evnnHOO+/9nnnnG/X/Tpk1d9+EjjzzSdQeuWbOme37XXXf11q5dW2A6fldiLQcNJejZs6frftulS5fI8ILXXnutwOeWL1/u7bXXXu49+g59p35rhw4d3HNbbrmlN2HChFANeyxsW0lmmJi6Ljds2NAt1wULFuR7/3777ReZ77vvvjvfaxpKpudPOOGEfM9/9913Xr169dxrLVq0cOvuwAMP9KpWreqeO+CAA7x169al9Ds11EWfveWWW/I9/9xzz0Xmr1u3bvleW716tetSqyGH0dtZUbbLQw45xL32yCOPBM7f2LFj3ett2rRxyzPWTjvt5F7/6quv0jYUtbCu+dHOP/98917tC9HOPvts9/yhhx4a+Lkvv/zSva519++//3q5prBld9JJJ7nXt99++8hzb7/9dmRb1PLQ+t9nn33ckG09N3To0KSHVQUdP9SdW++rUaOGe+2II45w/+//qRt6vM/65syZ43Xq1Mm9Xr16dbcfaj732GMPd9yP/cx1113n3tuqVSu3Leu92i4qV67snj/88MMDt8uiDHtUF3UN3/KHQ2jb0W+sX79+ZB5ih6L4v1XHA71f86X5HDx4cL51k2i4oYaQ6LXy5ct7PXr0cOe9nXfe2R23tO4ef/zxYp9LtG60jho1ahQ5VkWvO61bn/Z1/Q4dd/SbtIz9daahLK+88kqB+df3+8eJ6OnedNNNhQ7P0vrzywto+tpmtZ7bt28f2U7efffduOftK6+80v3eHXfc0X1OsYO/fO68804vlxV2nvP3T60HLdu6det6hx12mBvKuMUWW0T29X/++Sfw8z///LN7j5ZlsrQ+FFPocyNGjCjS7/KX/1lnneW26913393tE9quZ82a5d6jc45e889xen233XaLfPbqq6/ON81NmzZFhl1p/zzooIPcZ3r37h3ZFqKHWE2ePNm9z19G2o61b+k8qHOi9pVUPPnkk247059+h/Y3rRcNjdJ+GjQ8qij7drLHyIULF7p9StuF3qt9MHrfGz16dKFD6YpzzNM09T2VKlVy++zAgQMj+6zik2+++abAb9P607JQyQTF9loeWn7+Or/rrruyNuyxqPHd6aefHlnHfty17bbbuvXrr//YeS/O+VWIrYrmp59+imxr48ePL9I0Uj1uRa8vHSO0n2m9a3/RsUtxvF7TcT2WrkH8Zan36ru0HbZr1849H72PF/Va0d+/dFzq3r272/50bNX2r+8UHTP1nhdeeCFwmUyaNMm9rvhiw4YNkedPPPFE93zHjh3dfGuampa//LSNRXv00Ufddq/XNB/R23/0th4Ut2oopJalpv3nn38Gzuc999wTiRtjh+83b97cvda4cWO33+vayI+XdIxMpVzD33//HRn2G7u8Y69BVVLA9/vvv0eOW8p3aFloP6tVq5Z7TsfKJUuWpLy8RMdWf7nr92i70HnZj6WvueYaL51IfkXRgtZCfueddwJf1477xhtvFDjBaHz8ZZdd5j6rE2Y8d9xxh3uPdrh0nSR1sPQ3GD9QmTp1qvf1118XeK82yv3339+999Zbby3wun9wVAAbGxj4J6WgIPWoo45yr6muWWwSR8G9XtPBMHrcemlPfhW2rcSKFwzoZKHnlUzyKQGkg6QCFK1bHeSS+Yy/nZx22mn5amfogLX11lu714YNG5bS71RdAX1OB8Kgk54uoBUoK9At7DNF2S4//PBD93y8CwH/wBrvIkgJYb2uJEWskgjQnn322cj2F33SVe0//2JaJ6J4ySPtW6nw96tU/1IN3BMtO9WNUnCu15U0EP1GXRjpueuvvz5fQuj777+PXCjHJjlTSX4lW5Mk3md18eonS7U9xl6wq8Zi7P6uCxJdxMaaN29eJNHz0ksvFXjdX+5FqQOiIHDRokWR51esWOGCQb2mIDfot+qvWbNmcetCxkt+aX3oeTWw6FwTbdy4cS7oUXDy66+/puVckkzNLwXU0QFW9PPanxRAKwGfam2qeMkvv7aULrijkxfahv3fot8Zu73426Euvt96663A+dE+ETuvuaSw85z/O/SnpF50DRFdWPoXXbrgDqJlqGWn96hxMRk6n/nf+dlnnxXpd/mf14Ve0DlJF0valnRue/rpp/O9NmbMmEhA/sEHH+TbH/wLTsWKsXSci95v/XOojoextE1oeqlQIkjTi04U+xSXxV7kFGXfLsoxMpmaX/GO2cU95mm60cc8xaBqZPHnP5bWbVCdHjWeaVvRvhx74VoSya+ixne6XvETDUr6RFPMFO/8X5zzqxBbFS228hu2td9F/6ZkFeW4Fb2+9Hf55Zfnu1ZTfOMnPKMbkVUn1I8rgurqKt7Xe4p7rejvX/71RtB3Kcmi15XASpQoja0F9+mnn7r9J9a0adPc79Jnvv3223yvJXPtGi/eOProo93z0Y1u0XTu0OvR8YLOBX4DghqhovMPirX968HYzhCFOeyww9znlCwOauDw8wrRx08dh/WcEl4rV66MPK/zgN9IELtfJbO83nvvPbfNKsaKPe9pm/bXhdZXupD8iuLv4H/88UeRPq8eYNpgggKf6KRAbE+ZopwkleVVkOHvFPru6I0xHm3Ier8y6LH8A4yyz0EnX//iVTtG9AFOG62+P97vVkIwdocubvIr3l/Q9JI5Qem3ZXJbiRcM+Ce7448/PvLcxx9/7J5TIKD1pCDUPxHqYkHZfi3z6BO736tK6yGoV596Q+h1TUsBarLi9eJSckPLWgVHNd1Ro0YV2lusqNulkoBBQb0ulnSiV+tzvBa8hx9+OG5R5UsvvdS1OOkxUwGaDur+NhZ7sldPED1/22235XteJzS/JeWjjz5Kad781vZU/+KdjFPZnrVdqYi1WuP0mlqWlRyK7iGl3xzk9ttvjwQ+2Up+vf7665GWNV1cFZd/vFfvgVj+NpEsBZE6t2i/j71QFV2Q+T0Aoi9yogP22GC4sKSTLnT9Xs1ar0GUrA4KKotyLok3H6nwA8HYC/DiJL/8c2zQb9HxWAG5Xr/hhhsCt8MLLrgg8PvU2lycBE5JKOw8F538CiocrcBV26y23XjJLb/VXhfqyVAy1f9OXaQUhf/5a6+9NvB1v0dAbAu8Tz3GYht4lOSO7oFfGD8uitfyniqdB5ONY4q6bxflGFnU5Fc6jnlvvvlmgc/p4lmvKa5JpcC+38AdW2S9JJJfRY3v1ONNz2vegyjeykTyi9iqaLHVzTff7OZtq6228oqiKMet6H1Uyy2op7oSrrHHS8V2fiIkGUW9VoxOfsU7V2p6Ov4F9arSPq5eSvq8ehony792iL3RS3GSX/61XWxvxeikrdZ9dOLTb3xTwjCIjsMaQaRroaBGwXjeeustN101FsZ26NF61mtq/IotlK/lHJRU9m9gE3uuT2Z5+Um12F77sedWdXZIl4rpH0hZOq1atcr9icb8JqJbf6pOgsZG6zN+XReNZda/f/vtN1dXJZY/XY1HLwrVTPBrl0TTmGzVVVFNA5/G/2ustMb1//XXX24s7v+Sne716dOnx/0ejQ2PpfG3rVu3djVA5s2bZ82bN3fPayyupqkaAhr7G0RjtfU+zct//vMfSweNKQ8qVKkaMPFo7LDGLwdJpdBzKttKYTS2XLQ9+fzaEvvtt5+tX7/e1bHR7Zs1Ln3y5Mlu+9FYb9Wi8Pn1FnTrZ62rWKqFoDHxGsOtcdWaVjJUn2a33XZzNWu++OILV/NGNW5Ud+ykk06KzL/m2S+86c+//1q0omyXqqd36qmn2n333edqjPj8Wl7aJ1RvIEiife6mm25yf5kUXfMpdoy7lp/WhepP/Pe//408//LLL9uKFSvcWHnVkkqFtv94NYsywb/1eCwdC7R+VOsmevs87rjj4tYg0DJQDRzVTlStipL23nvvucejjjrK1d5JlmphqDCt9lPVwtP/a3vWOizsWJss1evRtqS6N0E1AlWvQsc31bPQvqp9NtYRRxyR0nfqWK91ofOLX8Mtll8rRvtzcc8lqdB8qaaL6n0sW7YsUtNvypQpkWUeXaeiODV2VHsw3rarfVrHH9W20nJXLYxkloGobonmX8sgF6Vynot3PtZ5SrGQakZqG9a+le64qDj69+8f+Lx/vIqu2RJ7vNL56PPPP3fnNNWm1L6pR9VRad++faSOVjyqK6W4SLVurrnmGneTENVxKipNT/OtW9KrELiWu+pypnPfLuoxMhvHPNX7DKqvoxjQj4VUYyc2JtRzOrb8/PPP7j2qLSc6N6XreJ6qosR3Oib66051dIJoPQbVSSwuYqvsKMpxK5quz4JqMelcJdHnKtW9VpynY5jqlGlb8us7BSnutaJqie2xxx6Bn9P0dCxXHU39XXbZZZHXtC+rfpaOj6pTFWvlypX27rvvumPiokWL3DWX6Pok3fu7apmqTpmmqbqD0fXa/Dp0On7r2BU9/xLvxgY6DquemZab9mXVNEzGQQcd5OJsnQdef/11VyNOtI787Ty6Jp6/bemYGn3tGXsDG+VHVLs13jEnlpa5atTpWjNerFRYnFkUJL/+RwG0L96OqUBQhQBVCC6R5cuXBz7v340oUQHYRHRCa9u2rfu3X8RQhWi1MUbvLDpJq2CpfyGQyjxKUCH16PmPLiro32ll5MiR7i+R2AJ+xXHYYYflK5iajEsvvbRAYc9MbSvJ0rJW8VKtMx0QVVRYySMtaxW71IW0ggg9p/UfL7Hkn5TinXx0QtNr2vaiT2AKDHTwiRV9ktd3KbjUdyv5FZ2cU5Cvi1f/OQWOEydOdBc0sQngom6XKtardae7meqEpAsKnaAeffRR9/pZZ50Vd3rF3eeKy1+2fkHOaIMHD7YLLrjAFd3UXXu0L0cXofQLBOcyXTz5CUkFUkpC6gR46KGH5ktIFrZ96r0qbqlimko4ZCP55RfyVVCXLP9uS0oGF+VYm6zClp/4d+YKSqboXJHqnRz9Y7uSP4Vth/GO7amcS5KlRIGCbf+iNFPLPHpZ6ngWdDfBwpZ7ppZBSUjlPJdou9RrSn5pv07HMbpBgwaRfyvZXNQbzki8m14Utr/561zrTuc87V96TsW3dSc6nZP0pyLwusDRhdyAAQNc3ObT+9SgpHOnYjgVi9axU8X5leTwGw5ESVIVfY+lY68u9OWBBx5w36O74elP60zT0EW+4tbo7bCo+3ZRjpHZOuYpToh390htc9reYvc9xRRKZPtJ30weW1JRlPhOsYf/++Jt5/GeLw3CGFv5xzbFQUHJqcIU5bhV1HOVji9K2GhZqiC+/rTP+deksQny4l4rFrat6sYhSnzp2iU6+eUnlYI6jrz11lvueS2Lktjf/SL4unb1b7QiimX8Oy/Gzqe/3HQM11+6rrErVKjg5uXGG290+4Wf/FKSS9+p9atzVqrHYyW/UmnQUyciJdzUESIosZ+pHALJr/+JvlBTy0BQoKsdSokvnfgViCiwUGuAH9Co5UnZXL8XS7xgMvZAnSwFOfEy+tGUAVeCQYGQ7lyj233r9ygQUNKgsA0sXmthotYXtfoqcEuksDuXhGlbSYWSS0oM6U4syqirxUrrTglNJbyUEddrutNHol5VRfHKK69EAtpEya/LL7/czYOfiNM2okSY/7oO5OrxqIscbf/+3R7TsV3qov3kk092dxPRnWW0HNTTUT0F1BIU1Cqcrn2uuPy7ZOqYEZ2gFt2pSctEFypafgoadDGiVjkt32T29aCAMLqlM1maPyUYU6WLr1xvDc0U3WFHSXhthwpY1INDjRPaphVYqIekLszjnQ9Kko4hqfKP7eoZoR4WicTrcZvKuSQZSoArcFRQrRZsHWeUKNXv0/FGPa90jMqFZZ6pZVBaz3Ppiot0EeQnytXSHa83QKb2i0R0hzNdRLz55psusaU/3SVUfzpv6dju9wbTeU3nVP0G9ahSq7b+fvjhB3dX1jPOOMPuv/9+996///47sIet+Mkv9c5QA5p6oeouy5qWvk//vvbaa90FpxqS0rVv57pU9zvFXephrmO37iaqXghKBmg96dii2EOv59KxpbhyMQFUlmMrvxem4mElEdTrMZf3GfUmV/yv452W7ZdffumukfWnuzLq+Kbev+m4VizsWK1GAyVfFHfp2KfrcTWOqEeUetSqQSGaEjRquFTSRdcj6qmkc4tiCy0HHUd1bEz3/q5tTw14upvk3Xff7X6XknDavrSdxjYs+MstXo+raGpsScXQoUNdvKT1pMYp3TXbTxZqeUWPJssU//dpuac6OqE4SH79j05wWtFq8VEWOCjQ8299+uKLLwZecPvdouPxs8uFbcDFoRbCSZMmuYy+DkCxJ4XC5jFV/pAVJWl0MVIWJLOtpEInD93qXUklXcjpYOAnt5QQUoJBPa/0XRoKoGSRhkjEDgGIbiWIl2GPfq8E3Yo4lrrU6kJI3YKVede86ATmD1fxk1+afz8giU3OFXe7PPPMM23EiBEuANUFrr+tJer1VVL7XDxqzfGPGfG6IqtbsQI0XRzp1vFajjrZ6v1FGQ6mLtzxLpIS0fZUlORXsrTNaRuIt33qAlgXs/57s8Fv9dR8JkP7ohJfClCDbhmdzmNtMvu3/1q6lp+//Wk/z5UEp78/qefXKaecUuD1dJ/f/GWp44hagIOO9ele7qXxPOefW4L45xgF1uk4RuvCRIkJHefU0q8eHummdamLZa3b7bbbLu4610WVEnHR9DvUWKM//3iiiww1jOoYG3t8ViOq38tLQ9U0BEVDX9STSxfwGiqjHuvJXITpvKrhvv6QX22zSqTpYkuJG/W81jot6r6d6jGyNB3zNCROy1gJTF0MZ/rYkoqixHdat4ofNXpADZxqbIyVTPyXi8IaW+m6Uj1rtC71Xakmv4pz3CoqJRqjeyXNnTvX7UMajqz4XEPgSuJa0e9VdeWVV7p1reTXs88+646papCILY2ihJMSXzomKtldUvu7ElRqtFOpGzXmKenmH4N1noil5abjrbbneMP0i6pNmzZu+1RvL21vWm/qWBA0L5k6Hvvbhdaf4uiSaiwsnU2SGeIfaNRVNoh/cRaUXX3//fcDh49FUw0BiVdjIR38eVQSJTbBIDoYpJPGDYsy/7k6fCMb20oqdCDUDq+kkrYjf0ihT4kkHcDV21AXImodiM3I+8M5lZgNWg9KOKlLvLqyprr9ad4UgCspp95XS5cuzTd/6gGmA5daD+L1TCvudqmgW71sND5dLUpq2dG0VOsi2/tcPOotp/lVslJDKYLoxKNhr7pIUTDnB1fRY+1ToZYrv4ZaKn/+eP5M8bfPeMGjnzzSsijuhYzfE9evA5Usvz7MCy+8kHDYS+w2HW+oQDqPtWrV1H6oIcVqEY6l4cB+PR7tq+ng92zWMS7RUOV0KmzdJToHq5VXx6CiTDceJWz8oSJBSYLo+hjpWu6l8Tynhg39xdJ2owYRbbvahmPpnPLLL7+kfIy+5JJL3HFV+4IubAujXglFOV7FSwz5xyv1Ogs6n0VTS77mV7T/JqJp6QLH741V2PsLo4Slekrqwk89VdUrojj7dqrHyOLseyV9zEt0bFFM5V8UZkNR4jvtH/6wqueffz5wulqPqSrq+kynsMZWfu9lUYO4aiElonWgYZ2ZOG4VlZIZSrbHHr9K4lpRyS8dM7S+dbxLNOQx0f6u9RZvn0nH9u8nlrSe1ICqmmPqARZU18tfbn6yN91O+l/PYc2LEsVabkqU+0OFY7ctHXOD6nOqc4TWd+y5vrDlpes4JX3Vu9w/npcEkl9R/BOoWuiC+EX/7r333nzPq6v5aaedVuj0/WJtqRZaTIVqMKnbtgqjxx50lelWTYp0Ul0ndVVUtl+JiKCWJAVJGs+cjYK22dpWUqHAVAGLTtJqydYFV3TXVz+R5LeWBA151NhsXYQrIFBLePSBRq1IF154ofu3MvtFKaobOw/RyS+1dKuVSd2L9V1quVJB63RvlyriK37tE7VkF3YCT7TPaRizlnN0fYB0UOuHWu5vu+22yDJL1B3ZPxGqBoy6HqvFtm/fvhYm6gWhCzFdCKvGQHQvBp00r7/++sgyKC6/h0mqCRvVKdPxTPuQ9qfYOhAK2BSkxJ4P1IIXmxxQD0VdqKSL9m3Nk5abtvvoedPxVb2gNH9q7Qwqdl8UurDQUC19p1pHNYQrluqSaFhVdACeyXXnL3MtX78wrd9zUAXpo+tUBU23KI0V/lCX6667Lt9FuJaLtlsFfDqG+z19wiTZ85yWhYb9Rtft0rrQc3pNMUJQbwutZ71P54dUkt7aDtSjSXS+0wWjf4OJaEr2qP6PbpqSCp1rdG5RL6zYJLaGw+hmHhI9DEr7gc6BsbXo9Pvffvtt9+/o84B6dgUVU9YQRw19jH1/Irpg0fIIqomixJ8arHT+9feDou7bqR4ji3M8Luljnn9sUZIkelvSd2gIaqLejZlW1PjO3+7vueeeAsdoDbnSjZRSlcz6JLYqXjJCCXAdRxRna3vUvhhN+4T2TW33SlgU57hVVIrbFOOo91QsxfQSvW5K4lpR26aWma6ldE5Q47f2m6D439/fVfrFL24vWtZ+A3u8umxK6Og47SfQUqXfr5hB61C92LU/a9kE9a7WcU7LUT1T1YgSdJ7TvPg1kFN1xBFHuHlR2RrVbYuXINYIJA1H1frW8VjnHJ86/ug5f7hk9Lk+meXlx/9KUvrbTuz2rmOVtuG0Sdt9I0NAt5zWIunRo0fg66+++qq7Vave07lzZ+/II490txKuVKmSe9xtt93i3qZdt1utX7++uzXz4sWLM3ZLZDn33HMjtxzVLYN1+/cddtjBPXfFFVdEbhsbK97zhd2GXreZ3Xfffd1rlStXdrdPHjhwoDdgwAD3bz2n13755ZeUbn+aaB50S+Rk+b/rgAMOSHhL4vfffz9t20os/1bCsbe+jb2Ntv5OOOGEfK/p1sO6Ha3/evRtvaPp1sP16tWLLNdBgwa5Wwf7twPX74+9pW2ypk+fHvn+atWqFbjd9vnnnx95/eSTT07rdhmtW7du7j3a53TL8kT++ecfd/tf3Vo5+tbBya6Twm7H3ahRo8i2c+yxx7pbPrdv3z5yjNDtlV988cVCpzd//nyvQoUKkd9/3nnnebmuKMtOt1b2t8WOHTu69a/jhtZR0Haf6JbRiY4f9913n3utZs2a7pbfuv23/qZNm1boZ2fNmuVuQ+3f0nn//fd387nnnnt6derUKfCZvn37Ro57eq/OCfpt2gYuv/zyuN+TzLYea9GiRV6XLl3c5zQvhx12mNe/f//IbbxbtWpV4PbzyR5n4x3bRbf69ud32223db9Zv3Ovvfby6tat657X7bhT+X3xvu/tt9+OLE/d2nvo0KFu3fnHvD/++CPynU2bNnW3vtZ+p+XRuHFj9/6g84OOezoO6DUdQ4YMGeKme+uttxbYpmPPszr+av/Wa9pWtc1qm/C3Ex0Px4wZE/e8HbtOCvu+XFLYec7fP7UOWrdu7dZNv3793H7nn4vatWvnLViwIPDzd9xxh3vPxRdfXKT5e/zxx70aNWq4aejYov1U60bzsM0220S2Q22vqe5/utW9zlV6n85TRx11lNezZ8/I8f3qq6/O9/4777zTPV+7dm23b+j9mg9/O9A2+uOPP0be7+/L2m8POeQQ7+ijj3bHEG1Pel4xZdB5K8i///4bObdqujouaDnsuuuukfm96qqr0rJvp3qM/Omnn9x86a93797uOK9974033ih0X8nUMS/o+7QM/ecVb+m7dHxp2LChV6tWrUj8Ensu0jFMz+uYlgr/u7Tf7LzzznH/xo8fX6z47pRTTnGvK8bQetW62m677dz/+3Hbfvvtl7bzqxBbFY+uF88666zI79X2eOCBB7pjysEHH+zOdf46vf/++4t13ErmXBQUh40ePTpy/tP0ddzQvukfG3QOf/fdd4t9rZjq/jVq1KjI+o533BMdW3fcccfItqzlqnnRNq/ri0suuSTu9+p36rXmzZu7/cnfBxItr1innXZavvkcO3Zs3Pf+/PPP3tZbb+3ep+OyjrVarzpGderUya1b7TNFdcYZZ0TmQ79d105Bfv/998hxS8dFLQedM3TO87e3JUuWpLy85O67745cC7Rt29atD/1GHZv0XXpe6yRdSH7F8BNYU6dODXz9s88+czuvElk68eskcsMNN7iTTqILiNdeey3uBV66k18K1keOHOl2bO3UChp23313d1BIFPgV9YJFNm3a5D3//PPuZKydUDuQDthaPvrNOlDqgJ7N5Fdhfwpe07mtRCssGPj4448j8/Hcc88VeF0nBz+wThQQz5kzxzvzzDNdQKUTiYI2BcAKYJMNpOPRgSsoUJJ33nknMv/xgpKibpfR/BOSDqCF8S+srrnmmsDXixugRf9pWeuYsP3227sLa63DNWvWJD1NXfz405o0aZKX64q67LSv6DPNmjVzxwidyPfee+/INpCO5JeORTfddJO7mPMvDKKPW4Ude1asWOHdcsstLhjT/lOlShX3XgXgsfOpY9ptt93mGkN0PtDFiS4GP/jgg4Tfk8y2HmTVqlXut3Xt2tV9n36fLvKHDRsWGHSkI/klSj7p4lzT0fLQctHFiIKvxx57rMB3F+dc8uijj7ogSr/Pn070uU+/SfPSokWLyLpRIPn3339H9s2g88PkyZPdOtSFk39xEB3YFnYBoPObnxTQtqvj4fHHH5/voi9sya/CznPR+6cC5lNPPdXt2zoeavmcc845CRv7lNjQuoi3jJKxcOFC7/rrr/f22GMPt24VQOv8othDF/7jxo0r8v73zTffuMB9q622ctNVTKOgXPt3rN9++81dWCo+1LapfXOLLbZw54RLL73Umzt3boFE7+mnn+6SsZpvLTMtO21jTz31VL54qTA6tz/00EPuvKjku86tuiht06aNS+Aovogn1X071WOkKP7TBbLe61+ER++jifaVTBzz4n2ftiVdCGq56TcpYX7MMcd4M2bMiHsuKm7yq7C/6GNkUeI7xV3+MVXLTscvnaN0LfP0008HxlPFOb8KsVV6TJkyxSVddZzUetMxSMcUJUW1/f/666/FPm4VNfmlxuebb77ZXfMpCa19U9cnSsZoG413Xkz1WjHV/UsN836SWMcaNZjFo+OYlqMSdtqWlWTRce+HH35I+L06p+lcp+O85j/2fJJM8kvJbP9zSmxpP01EiUM11ml/92MQJUF1DFYjxldffeUV1fjx4yPzooR2Ivrt6qyhY7CWmda7zmHaFlavXh33M4mWV3SMpnO2Gsz8aetYp8T+Pffc482bN89Ll3L6T/r6kZV+6gKpLsbqWqwC2+miAq3vvPOOG/ajYuEo/TK1rSCYuiOr/o6Kt6pLsl/PIogOa7qjjLryqqu87moFAMiN85zurqebqWjYnQr/AihZGhaoukjarzNx4wgAyEUkvwJobKvqeOiuGem4S5xuZd2jRw9XjM8vwIdwSPe2gvhUAFS1N5T0ijce36fx8brDi4r9qrYJACB3znMHH3ywu0GK6rKoMDWA9FNtLhVrj75Jkm40MXLkSFenR3eEVANh48aNszqfAFBSSH7FKeSnFkkVa03HLVlVLFx37VDxVXqghEu6txXkp4LAKm6qYom6E4gOVyrgm6jArXqIbbvttq4wo24RrLuoAABy4zynAuu645gKUesOwgAyQ43uulOcCo7rphIqKq4bf6jguG6CoELZQXfDA4CwIvkFIGfpzpC645juFqK7B6knl4bJAAAAID7dfVMJLg0z1l3ZdGe5hg0bWs+ePe28886zXXbZJduzCAAliuQXAAAAAAAAQqt8tmcAAAAAAAAAyBSSXwAAAAAAAAgtkl9ABooDq5Do2Wefne/5J5980sqVK+cKkKZChUn1Od2xJ100LU1T0wZKu2xvz8uWLbMtt9zSdt55Z3dTBgBAbsdkMmbMGFdL9JBDDrEmTZq484j+/vzzz6zMK8Jfx1bb11577ZW1eXj22WfdPDzwwANZmwcgm0h+AWmmAEt3GLzyyiuzPStAqadksQI1JY9zVZ06deyyyy5zd/V9+umnsz07AIAkYrKjjjrKrrnmGnv77bftr7/+ysr8ITz85Gku0zbfuXNntz8sWbIk27MDlDiSX0AavfLKK/bll1/aGWec4e6okw66PfUvv/xiH3/8cVqmB4SN9g3tI9pXsuWss86yBg0auCTYunXrsjYfAIDkYrLDDz/cbrzxRnvvvffsn3/+yco8ouzo0aOHi1Wy2UhWvnx5Gz58uEt8XX/99VmbDyBbKmbtm4EQuvPOO93jiSeemLZpVqpUyTp27Ji26QFh06ZNm2zPglWtWtW1qN5999324osv2pAhQ7I9SwBQphUWkz3++OMlPEcoy6pXr54T8fyhhx7qGutGjhxp1157rdWsWTPbswSUGHp+AWmsK/HVV1/ZLrvsYh06dEj43lWrVrkeIm3btrUqVarYVlttZccdd5zNmzcv5ZpfP//8sx1xxBFWv359d2JVd+a77rrLNm/enFQtpE8++cT2339/22KLLdzQgB122KHQVim1ph544IHu5Fm5cmXX4+aYY46xqVOnJpz/TZs22R133GHdunVzJ9tc7B6+evVqt/x23313t0y0flq2bOlqgjz//POB77/55pvdcqtVq5ZbB9tuu61dccUV9u+//yZcHqoP9cgjj9iOO+5oNWrUcMPntC6+/vrrwHmbMWOGDR061Fq1auXmS8tQ83bwwQfbE088EfiZX3/91U499VSXIFKCRt+x5557uroPQVSLQvOn2hSff/65+91az2ot1NDDwYMHu9f1m+PREBK9R+vZt2HDBvedRx99tAv+ateu7bY37SvnnHOOzZ8/P3A5PfXUU+7/TzjhhMiQAv2pTosvdjvXdtasWTP33DfffBN3Pv/73/+695x//vmBvcnUK6Bx48ZuG1evgX79+sVdN+LX87v//vvjvgcAkFsxWWml87t6ten3KfbQebVTp07uOcWGsaZNm+bOpYobFEPUq1fP9t13X3vppZcCp6/zrH++XbhwoZ155pnWvHlzd07Uo4aULl26NPCzL7/8svXu3dvVw1Qjrh41byeffLJNmjQp47Gl5ksxhuq9BcXWvv79+7v3q+HKN3v2bLvllltsn332sRYtWrhlVbduXRcXPvzwwy6+DlpOvuhYJTo2Car59f7777vnttlmm7jzuHHjRnedoPf99NNP+V5bs2aNjRgxwm3nmkfFedoeLr74Ylu8eHHg9LQ+1Fi3fPlye+aZZ+J+LxBKHoC0uOqqq1Tp2rviiisCX3/iiSfc64cddpi3/fbbe3Xr1vUOOeQQr2/fvl7Dhg3day1btvSWLl2a73MzZ86MvBbr008/9apVq+Zeb9OmjXfkkUd6++23n1e5cmVv0KBB7jN6TdOI5j9/5ZVXeuXKlfN23HFH99lddtnFPa+/O++8s8D3bdiwwRs4cKB7vUqVKt5uu+3mDRgwwOvSpYt7TvPy7rvvBs5/ixYtvEMPPdTN27777usNHjzYLYdcMmfOHK9Tp05ufqtXr+6WpZbLHnvs4dWpU6fAOli8eLHXtWtX9/7atWu733fEEUd49evXd8+1atWqwLKPXp/HHXecV6lSJW+fffZxy7V9+/aRZfvNN9/k+9zkyZPdd+j1Dh06eIcffrhb9rvuuqtXs2ZNtw5ivfTSS17VqlXdZzp27Oj169fPfVeNGjXccyeccEKBz/Tq1cu9dsYZZ3jly5d3y0PLYP/99/eef/557/33349MLx7Nm95zzz33RJ6bO3eue07LUduZ5r1Pnz5ekyZN3PMNGjTwZsyYEXn/woUL3fLRdq3Xe/bs6f7f/xs9enTkvUHb+WWXXeaeO/XUUwPnUdtyo0aN3HsmTZqU77ULL7zQPa/f36NHDzevO++8s9tXKlSo4D3++ONxf7t+hz47f/78uO8BAGQ3Jgvixz86X+W65557zsUKfnyl2EPneMUCOlcNHz483/vffvvtSDygGELndcUDOqfpuaFDhxb4Dk3Df61Zs2bunKnzu87dOpfrte7du3vr16/P97lrrrnGvVaxYkVvzz33dPGePrPddtu5eYuNLzMVW+rfes9NN90UuAwXLVrkPqc//dt33XXXRWI4TVPLSrGR3qfntQw2b94ceb/iEcUl/vYTHavoT/GMfPLJJ+51Tcu3adMmt2z1/Ndffx04n2+++aZ7fYcddsj3/Lx587zOnTu71+rVq+f17t3bbQN+TLT11lt7s2bNCpymtge9R7EdUJaQ/ALSZPfdd3cnknfeeSdh8kt/BxxwgLds2bLIa0uWLIkkUW688cakkl+rV6/2mjZt6l7TxbpOoL4pU6ZELuwTJb+UeHnrrbcC51OBjb4j2rBhw9xrSgT88ccf+V57+eWXXRC1xRZbeP/++2+B+defTvDTp09PYmnGX3ap/MUGf4lo+e20006RYOCff/7J9/qaNWsKrFslGP3lER04rVixwjvooIPcawriokUvD62H6OWxceNGF2QGBSRKVOn566+/vsC8az2NGzcu33NK6CiIVLD76quv5ntNwZAfMD311FOByS/93X///YHLScFmvEBNQZ62q9hgcvny5d4bb7zhrVu3Lt/7FTT7iSoFx7H8gFLbQDxBya9ff/3VPacks9ZdLM2LXlfiN9ojjzzinm/btq33008/5XtNy7hWrVrut2n6QRSE6/PPPPNM3PkFAGQ3JstU8iv6HJrKX2yclsgPP/zgzrNKJKmRKTr+88/xeo/v77//jiSrFENEJ26+//57F7fpNZ3/gpJf+jv++OO9tWvX5mss9GNQNYz59B4lq9QoN23atALzrnn75ZdfSiS2/PDDDxM21t19993udSUOo3333XeuwTGWkk1+Qk6Ni7H8+YknKPkll19+ecLGOiW09Pq9994beU7rUI2Cev7EE090MVZ0MtFvxNt7770Dp6nGW20/auiNjcuAMCP5BaSJ35sm9sQdm8DR+4J6hYwaNcq9rpa4ZJJfTz/9dOT52FY3ue+++wpNfl1wwQWB86pAQa9/9tln+U6UCmiUTPnzzz8DP6feQrEn6OgARfNcFJ9//nmBlrRk/qJ7BxXm9ddfd/PYuHFjl7wqzOzZs13PIAUPsUkS0TLyW1m//PLLwOWh1rxYf/31V6T1M3q9KjGk5ydMmJDU7/ETc7fffnvg6wrugpI/fuAeux1GU49BveeUU04p8Npdd93lXuvfv7+XCvUA0/KMDuCKk/wS9diLDcx96oGp17Sf+HQB4fdEi75wiHbrrbdGEs5B/ETe+eefX+hvBgBkJybLVPJLvYyKEq/4vYOS4Z+/zj777KTe7/dkij3f+xQn6PV27doFJr+UXFq1alWBz918880Feo2p4VDPJduzP5OxpRJEfnzw1VdfFXjdb3RWL6hk+b3f1TMtXcmv3377LdLoHNtYp+WpRKdiQi0rn3rC6TP6DUp2xVI8o552ek9QIk8U7+r1oBgWCCsK3gNpoBpe+hPVNUhkp512cnWEYvnj/RPVJog2btw49zhgwAA3fj+WaivpDnSJqJ5TEM2LakNEz4tqg6m2gOpDxLurnuoYPPDAA67ORtB3qzZZUajOgv4ySXd7EtVBSKb452effebqPqjW1/bbb1/gdS2jAw44wN544w237Hbbbbd8r1esWNHVtoilug6qNaZ6YarXoP/37xI0ZswYO/30092t2Xv16uVqOwTRfL377rvu34MGDYq7Hep3qi7K2rVrC0xLdTAS1bbSXYJU2F310VRXw+fXHlNtsiCqV6F6WjNnznT7jF87QzUt9O/ffvstX62w4lBtE9Ut82uV+VS75J133nF1PLS+fVoWqj2m+miqwxbEr9WhbTyIv/8vWLAgLb8BAJC5mCzdLr300oxOX/WtPvzwQ/fvU045JanPqNaUqLZsEN0QQDUwVVdU58AmTZrke11xn2qKJRO3ql6X6nCprteFF17opq1aX/FkMrZUjSz9ZhV1Vxyw6667Rl6bOHGi+1M8HhSL6a7NH3zwgX3//ffuTqD6f+W3VqxY4V6fPn26pYtiDtViVVw5evTofPHKc88952qmDhw40NVo8ymG8X+74slYqtOqaar2m5bbdtttV+A92jf++usv4hWUKSS/gDRYtmxZ5N8qep6IimcGUaFSUSIiGX/++ad7jFcIX4UvVdw8et6KMy9//PGHe1TiorBC9UouxFLB8KDgKVeowKkkeyceP9hT8fnC7kIYlNBUwBWUtPSXv5Jf0cv/oosusi+++MI++ugjF6jps126dHHBzZFHHmndu3ePvFdJMxUyFRWlLYzeHxt0xtuupHXr1i75poBagZqfQFLySMktBc4q3B9NFyLHHnuse38i/nyng4JFFdPXMtP+oiL4osL7CiaVGFSiMXYb//3334u0jUfvO0E3OwAA5FZMVtrofO0n9pIt5F9YvKJ4UYmVJUuWuHNlbPIr1bhVN01SA5qK0OtP0955551tv/32c3GAbtBUUrGlGsGuu+66Ao11fkOd7sysovjRdKMcxQdz5swpkVjFbzBU8kvzFZ388udTvyOav9yuvPJK95cI8Qrw/0h+AWmgwMGnViH/hBJErTHplChYKCyQSGVe/B46ukNlz549E743KIEU3TsoVUr6PPbYYyl/7rDDDnN/uSjV7UDBnVp71QqpXmpqydPfDz/84IJL3d3Jv8tg9J2I4rX0RlMPqFTXlwI1Jb/UmuonvxIFk7q7qRJf2jZ0p0gl6xQA645Oop5xupNi3siB9NAdNJUA0+3sFYwPGzbMPa95Dgom/eWm3nbqtZdIdPAedNEVnVQDAORmTJZuOr+p53yqbr/99rjnldIWr+yxxx7uDofqnaRRCopVdFdD9UgfPny4iwXU06skYks15O299942duzYSGOdGr/8u3fHxgG6g7fiRvWG0mvqba950zakuEZ32FTSMZ2xij+KQ3fPVBLQb6ybMGGC60GnxsnYBkV/uWlUhN/QGo/uQB6EeAVlEckvIA2UmNCFtlrj1CpXEoGW31PHv4Vy0Ekt3i2oi8LvQaSTvp88KCkaCvfUU0+l/DkFPckmv/yWzWSDVn/5+61vQfzX4nXlLwoljfxeXhoq+Prrr7tkk4YEqKVVQZ4CaAWEGkqQqYBaXe01/ECB2ty5c61Ro0Zxg0nxb6Wu1tegYaIabpEJmhclv7TNKvnlB5MKLNUKHbSNayhAUbdx/9biWh4AgLIRk/nUOOWXpUjF1VdfndS5Wucn/T4laTT0Lmg4WyzFIIpt4sUrihfV68t/bzooBlFM4pdQUO+jK664wh555BHXeOb3ti+J2FJxgJJfaqBT8uutt96yRYsWuUa32N5z6n2lxJdKWih2KKlYRetUjXUjR4508e7ll18eWR5qxIxNQPrLrW/fvm7IalEQr6AsSm8XFKAM04lSpk6dWiLfp+Fu8vLLL7skSCw/EZEuaqVTLx319lH9g5KkGlP/u0FHSn8KJpPl13x44YUXIkMKClv+CkZUM0JD/WKpjoJfR0wJqUxQnQcFln4vJc2LqHXST+z4SadMBGoaFqDWR/WqUjCpQEott+3bty/wfj+wbtmyZYHX1CKsQDSI3zMsaBtPhlpFNT8KWL/88stI77SgYNLvjaZ9eMqUKUX6PtXXkHg1wwAA4YvJfIqRihKvJCo1EC36/P7oo48m9Rm/VmW8RkQ/ydOuXbu0NtZFUy2wW2+91f1bwwn9oXYlEVuqsU5lQJQAU2NdvKGE0bFKvKGeKpsQj1/Koqjxil8rVetJNcb8OF4xcKyDDjoocg1QlF5oitf+/vtvF8v5tduAsoDkF5AmfoJDQ7dKgrpIq26Uen6phSh6qJta+FTgM53UMqQu2UoMqVD+5MmTC7xHJ+s333yzSF3+s+3QQw91hdZV7FXL1m8R86mmhV9E3g+M9D4FHaeeemq+92sZqRCtPqOWxdhi90Whnl1BBVYVvGjoY2xiSUMLFFCqVpgCqejtIzpR89prrxV5nvxATa2TfvAcFEyKH1zde++9+Z7XbzrttNPifodfp6uoyajoeXrooYcSBpMKXLXctE779evnhtsGFRtWAK2aIEH8/X+fffYp8vwCAEpXTFaSFPOp8eu+++5zsUFs8kO9qsaPHx/5/5NPPtn1flPP5xtvvDHf+1WrUzewEcULxaXvVpmKoJpYaiTzh9n5vfFKIrZULzTVRlUcdMstt7iGSb8BL16sol7tsYlT9VpT7/VMxSt+TzQ11l1yySUurlQDnpKSsdTjSw123333nYtxgup6KcGouCcoGefftEfTj1d/FgilbN9uEgiLCRMmuFsG9+jRI/D1J554wr2uW1oH8W/brNsyJ/O8fPzxx+720Hq9bdu23pFHHuntv//+XuXKld1tmFu0aOFemzdvXr7P+bd+1rSDaB71uuY5mm6nfNRRR7nXypcv73Xr1s074ogjvEGDBnk9e/aM3Fpct2BOZv5zzaxZs7wOHTq4+a1evbpbloMHD/b23HNPdwvq2N+waNEir0uXLpFbVOv24/379/caNGjgnmvVqlWBZZzM8ghaP/73aJqHHHKId/TRR7v50y3C9fw+++xT4HbXL730kvsd/q3K9X597qCDDnL/r+e17qLpFtx6XrfkTsY222wTub231v+KFSsC3/fqq6965cqVc+/r3Lmz21Y1z7qFtx532223wO/VLbi1remvd+/e3gknnOCdeOKJ3htvvJFweUXT9l+hQoXIfGp9JnLRRRdF3rvtttt6ffv2dfO71157eXXr1nXPP/jggykfAwAAJSOZ4/G1117r7bzzzpE//7iv2MZ/7vTTT/dy0VNPPeXOn348odjj8MMP97p27erOtcOHD8/3/rfeeisSL3bs2NHFNvvuu69XsWJF95zOrbE0Db0WOy2fztd6XXGD78cff3TPad66d+/uDRw40P1pmep5zdtjjz1W4rHlN998E1m/+hsyZEjc9+qcr/collbcpPO/lpnm/fLLL4/7vf/973/da/Xr13e/WbGK/hQrxltesW6++eZ88/n444/Hfa9iG61vP/5SHKV59bcDP+5Zs2ZNgc+ec8457rUHHnggiaUHhAfJLyCN/Av4qVOnlkjyy08O9OvXz6tXr54LbDp16uTddttt3rp169yJW4FE7ImvqMkv35gxY9zJtWnTpi7AUUJASRCddJ9//nlv1apVSc9/rlHy5pZbbnFBW61atbwqVaq4eT/00EO9UaNGFXi/futNN93kAg0lmrQOtCyGDRvmLVmypMD7i5r8evvtt10QrqBQyTWtWyWwlJBRELx+/frAaWka559/vrfddtu54Ejzp+nrcwqyfvvtt2Ilv2699dZIkBZv2/Z99tlnLthWYKhlpXm64YYb3Laa6HtHjx7tAmCtDz+BFh2MF7Y9S58+fSLzGW+7jvbll1+6RKGmrW1A392+fXuX4FTgHrRu/WBS6wMAkLsxWXSsk+gvUaIi26ZMmeKSK2oU03lKjXCKAc866yz3WiwtB/1mxQ5+7Lb33nsHxjZFTX4tX77cu+uuu1xc2q5dO69mzZou9tD5UwmnH374Ie7vyXRsqcYsf70minEUTymOVkOdYhXF10qCffDBBwm/V7H2xRdf7BqjFaP53+XHJskkv+bPnx9JWiVqUPStXbvWe+ihh9x63HLLLV0ys2HDhi4mPfPMM733338/8PcpDqtdu3ah0wfCppz+k+3eZ0BYvPLKK24o3AUXXGAjRozI6ryoaGevXr2sc+fOrsA3gMzREFcVoNXwgZkzZwbeQRMAUDZjMiBXvPrqq65e7Pnnn+/uFg6UJSS/gDTT+HkVHv/9998zfgcVjfFfuXKltWrVqkAtJxX41C2ZdWLTCQ5A5tx222128cUXu/pnKqYPAChbMRmQ61T3rGvXrjZv3jxXW6xevXrZniWgRJH8AtJMxUN32mknO/30010x0kzS3XFU1LVTp07WunVrV9RTvU5U1FQnON0RaMyYMa4wKoDM0G3itf+1bdvWFcIvV65ctmcJAFDCMRmQ63S3ymOPPdbuv/9+O+OMM7I9O0CJI/kFlGK6M6Hu3DNu3DjXirNixQqrVauWbbvttnbUUUe5O/yQ+AIAAAAAlGUkvwAAAAAAABBa5bM9AwAAAAAAAECmkPwCAAAAAABAaJH8AgAAAAAAQGiR/AIAAAAAAEBokfwCAAAAAABAaJH8AgAAAAAAQGiR/AIAAAAAAEBokfwCAAAAAABAaFXM9gwgvs2bN9v8+fOtVq1aVq5cuWzPDgAAoeN5nq1YscKaNGli5cvTJlgaEB8BAJBZXgjjI5JfOUyBXfPmzbM9GwAAhN7cuXOtWbNm2Z4NJIH4CACAkjE3RPERya8cphZNf4OrXbt2tmcHAIDsmTjRrFcvs3HjzLp2Tdtkly9f7hIp/jkXuY/4CACADJo40Zb36mVqZgpTfETyK4f5XfkV2BHcAQDKtJo1//8xA+dEhs+VHsRHAACUQMxl4YqPwjF4EwAAAAAAAAhA8gsAAAAAAAChRfILAADkvipVzNq0yXsEAABAZlSpYtaqlYUNNb8AAEDu23Zbs99+y/ZcAAAAhD/mmjjRrE4dCxN6fgEAAAAAACC0SH4BAIDcN2mSWYMGeY8AAADIjEmTQjnskeQXAADIfRs3mi1alPcIAACAzNi40WzJEgsbkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAAAAAAAgtkl8AACD3tW9v9tVXeY8AAADIjPbtzT780MKmYrZnAAAAoFA1a5rtumu25wIAACD8MVePHhY29PwCAAC5788/zS64IO8RAAAAmfHnn2aXXWZhQ/ILAADkvn/+MbvzzrxHAAAAZMY//5g98ICFDckvAAAAAAAAhBbJLwAAAAAAAIQWyS8AAAAAAACEFskvAACQ++rXNzvjjLxHAAAAZEb9+mYnnWRhQ/ILAADkvhYtzO6/P+8RAAAAmdGihdmIERY2JL8AAEDuW73abMKEvEcAAABkxurVZhMnWtiQ/AIAALlv2jSzHXfMewQAAEBmTJtm1quXhQ3JLwAAAAAAAIQWyS8AAAAAAACEFskvAAAAAAAAhBbJLwAAkPvKlzerVSvvEQAAAJlRvrxZzZoWNhWzPQMAAACF6trVbPnybM8FAABA+GOuefPM6tSxMKH5FAAAAAAAAKFF8gsAAOS+qVPNtt027xEAAACZMXWqWY8eFjYkvwAAQO5buzYvGNMjAAAAMmPtWrPp0y1sSH4BAAAAAAAgtEh+AQAAAAAAILRIfgEAAAAAACC0SH4BAIDc17q12Rtv5D0CAAAgM1q3NnvhBQubitmeAQAAgELVrWt26KHZngsAAIDwx1x9+ljY0PMLAADkvr//NrvpprxHAAAAZMbff5uNGGFhQ/ILAADkvvnzzYYNy3sEAABAZsyfb3bttRY2JL8AAAAAAAAQWiS/AAAAAAAAEFokvwAAAAAAABBaJL8AAEDpuPNQ//55jwAAAMiMunXN+va1sKmY7RkAAAAoVOvWZi+/nO25AAAACH/M9fTTZnXqWJjQ8wsAAOS+9evN/vwz7xEAAACZsX692bx5FjYkvwAAQO77+Wez5s3zHgEAAJAZP/9s1qmThQ3JLwAAAAAAAIQWyS8AAAAAAACEFskvAAAAAAAAhBbJLwAAAAAAAIRWxWzPAAAAQKG6djVbu9asUqVszwkAAEC4Y65//jFr2NDChOQXAADIfeXLm1Wpku25AAAACLfy4Yy5GPYIAABy36+/mu21V94jAAAAMuPXX8369LGwIfkFAABy38qVZuPG5T0CAAAgM1auNPvySwsbkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAAAAAAAgtkl8AACD3tWhh9uijeY8AAADIjBYtzO65x8KG5BcAAMh99eubnXRS3iMAAAAyo359s+OOs7Ah+QUAAHLfokVmjz2W9wgAAIDMWLTI7KmnLGxIfgEAgNw3Z47ZySfnPQIAACAz5swxO+ccCxuSXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAIPfVrGnWq1feIwAAADKjZk2znj0tbCpmewYAAAAK1b692aefZnsuAAAAwh9zjRljVqeOhQk9vwAAQO7bvNls3bq8RwAAAGQ25goZkl8AACD3TZxoVrVq3iMAAAAyY+JEs4YNLWxIfgEAAAAAACC0SH4BAAAAAAAgtEh+AQAAAAAAILRIfgEAAAAAACC0SH4BAIDct912ZnPn5j0CAAAgM7bbzmzqVAubitmeAQAAgEJVrmzWrFm25wIAACD8MVfTphY29PwCAAC5748/zAYMyHsEAABAZvzxh9mQIRY2JL8AAEDuW7rU7JVX8h4BAACQGUuXmr3xhoUNyS8AAAAAAACEFskvAAAAAAAAhBbJLwAAAAAAAIQWyS8AAJD7mjQxu/HGvEcAAABkRpMmZlddZWFD8gsAAOS+rbYyu+yyvEcAAABkxlZbmV14oYUNyS8AAFA67jz05pvc7REAACCTli41GzPGwobkFwAAyH1//GHWt2/eIwAAADLjjz/MBg+2sCH5BQAAAAAAgNAi+QUAAAAAAIDQIvkFAAAAAACA0CL5BQAAcl/VqmadOuU9AgAAIDOqVjXr0MHCpmK2ZwAAAKBQSnxNmZLtuQAAAAh/zPXdd2Z16liY0PMLAAAAAAAAoUXyCwAA5L6JE81q1857BAAAQGZMnGjWtKmFDckvAACQ+zZvNluxIu8RAAAAmbF5s9nKlRY2JL8AAAAAAAAQWiS/AAAAAAAAEFokvwAAAAAAABBaJL8AAEDu69jRbPz4vEcAAABkRseOZuPGWdhUzPYMAAAAFKp6dbMddsj2XAAAAIQ/5ura1cKGnl8AACD3zZljduaZeY8AAADIjDlzzC680MKG5BcAAMh9ixaZPfBA3iMAAAAyY9Eis8ces7Ah+QUAAAAAAIDQIvkFIKGPPvrI+vTpY/Xr17dq1apZx44d7fLLL7eVK1emZfo//fSTVa5c2cqVK2dt27YNfM+sWbPc64n+Lr300sDPfvbZZ3bjjTfaEUccYVtvvXXk/V988UVa5h8AAKC0S2e8t2rVKnv++eftwgsvtL322stq166dMM6L5nmePffcc9a7d283L5UqVbK6detaz5497Z577rH169cn/PzMmTPt7LPPtvbt21v16tXdd+u3nHDCCfbHH3+k/FsAhAcF7wHEdeedd9oFF1zgApY99tjDGjVqZJ9//rlLJr366qsugaTApKgUwAwZMsQ2btyY1Ptr1Khh/fv3D3xtxx13DHz+nHPOcQk2AAAAZD7emzFjhh199NFFmpfBgwfbiy++aOXLl7fddtvNmjZtagsWLLAvv/zSvvrqK3vhhRds7NixLkEXS68NHTrU1q5da507d7ZDDjnE1qxZY7///rs9+eSTriG0devWRZovAKUfyS8AgX788UfXYlehQgV766237KCDDnLPr1692g499FD7+OOP7bTTTrNXXnmlyN9x7bXX2qRJk+yss86y++67r9D3K/BS8JKK/fbbz/r162c77LCD+1PL4ezZs4s8zwCypGFDs/PPz3sEAORsvFerVi3X00pxV7du3Wzp0qX2n//8p9DPjR492iW+6tSpY+PGjbMuXbpEXlOvrT333NO++eYbu/vuuwv0+Nd8HnPMMdawYUN76aWXXBIvdhRBlSpVkv4NQJnWsKHZGWfk1VoNEYY9Agh00003ua7nCl78QEjUhXzkyJGuRU6tgdOmTSvS9L///nu7+eabbcCAAa4lLlNuu+02Gz58uGv9U+shgFKqWTOzO+7IewQA5Gy816ZNG3v88cdd46YaHdVzPxnq0SXqNRad+BL12DpDF+Nm9vXXX+d7bdOmTXbSSSfZ5s2b3bzGJr5EpS8aN26c9G8AyrRmzXRwsLAh+QUgcDjiO++84/591FFHFXi9ZcuWLpjxW+lSpe7oxx13nG2xxRZJ9fgCAFPdGV3wpKneIACUdZmO91JVtWrVpN4XOwRTPdbUs2v33Xd3QyUBFNPKlWbffWdhQ/ILQAG//vqr6+4uO+20U+B7/OfVXT5VV155pf3yyy+ucKm6p6dSQFW9xdT9XsVM77rrLps+fXrK3w+gFPr1VzNd1OgRAJDz8V6q/J5nKpYfW69Vwx4ffPBBV5fs5JNPzvfa+++/7x41LFJ1ZNX7SzXMTj/9dLvllluKPEoBKLN+/VW1YyxsqPkFIPBOOaK766huQ5DmzZvne2+yVKz0jjvusL59+7qipqlYtGiRXXbZZfmeU3Cj7vEKiGrWrJnS9AAAAMqqTMZ7RbHPPvu4O0zecMMNkTqtfsF7Fd3XvLzxxhu2yy675Puc6sdKxYoVrUePHgUSdcOGDbPzzjvPbr/9dpc8A1A20fMLQAErVqxwj4lqNPiJpuXLlyc9XbUuHn/88a6QqZJVyVKBUrXyqWVv7ty5bjpTpkyx6667ztWkePbZZ13dMNWsAAAAQPbiveK4/vrrXVyn+E53nBw1apR98sknLsbr3bu3bbvttgU+s3jx4kj9sjlz5rjPL1y40MWMqv2qpJgaXjV6AEDZRfILQInRnXl0+2sNV0yl6Kje+8gjj9j+++9vzZo1c7e37tSpk11xxRWuOKruUPTBBx+41kAAAACUPhs2bLChQ4e6uzYedthhNnnyZFfyQsMzVTz/0Ucfte7du9vEiRPzfc5v/NTnNWRSIwJUF0wx43//+1/XWOonxzQ9AGUTwx6BMkrBgIYRxnryyScjXd8TBQgr/1d0unbt2kl936effuqK2/fp08eGDBli6aLu7bqT4+uvv+4KnipYAhBCFSuqynHeIwAgJ+O94rj11lvtiSeecLHiM888E3m+Xbt2duedd9qaNWvs4YcftnPPPdfGjRsXed3/HbqjoxpKY6n21yWXXOJ6un333Xe29957Z/y3AKVaxYpm9eqZLVliYUIECZRRr7zyis2ePTswGFLwIEuXLnWBQlAdCHUlF/+9hVFySi1z6o6+11575XtN3yPz5s2LvKbeYV27dk1q2ttss42b/p9//pnU+wGUQttvb7ZwYbbnAgBKlZKO94pD8yTxasLqjpRKfqn+17p161xZDGndurWNHz/ePQbR72rQoIEbCvnXX39l8BcAIYq5Zs40q1PHwoTkF1BG6ZbQ8XTo0MHVWlBtrR9++CGwhUzPiwqSpuLnn3+O+9ratWsjLXl+QiwZfq2HeMVaAQAAyqJsxXtFoQbSRL3MVDNWNm/e7OLERo0auf/fcccd7eWXXw7s4SabNm2KxJXcHAkou6j5BaCAypUr28EHH+z+rdoJsdSCqLs2Sr9+/ZKapnpyqedX0J8KmUqbNm0iz8X2DotHXfU13NEfAgkgpKZMMWvbNu8RAJCT8V5x6M6O8u233wa+/s0330QaO1XTy6ebHukujtOmTQscBaDSG6oHpvfstNNOGZt/IDSmTDFLcgROaULyC0Dc4vQKElR74b333os8r9bBE0880bWiKdjo2LFjvs+ploKei32+OFTs3u92H0233e7bt6/rwq7bdKtIKoCQWrfO7Pff8x4BAKGL9/r37+8eVd9LCatoKnJ/5ZVXun8PHDjQ3ezI17ZtW1ckf/369e7u4MuWLcuXwDv77LMj02/SpEna5hcIrXXr8oY9hgzDHgEEUvf2ESNG2AUXXOAKj/bq1csaNmzobjutZJO6yj/00EMFPqdgafr06WmdlwceeMBOO+00d3vr9u3bu5ZKJb4UCKnmw5ZbbmmvvfZavlZA32OPPeb+fH6th1NPPTUyTFJ3kxw9enRa5xkAAKCsxnvqKebHXMuXL3eP6pW1yy67RN5z0kknuT+fkltKen3//fduCKbu7NiqVSv3OfUGUyKuc+fOdssttxT4vnvvvdemTJniEnhKhul7VE5DvcVUtL9Lly6BvwNA2UHyC0Bc559/vgsyFBSphU9DDFu0aGGXXXaZ+yupGlvnnHOOvf/++zZp0iQXFCmIUs2G7bff3gVqZ5xxhgvUgvgBU6ypU6dG/t2yZcuMzj8AAEBZivd+/PHHAoX21WAZHZMdeOCB+V7X96iYvZJUquGlZNaECROsRo0aLhGmHmhnnnmmVatWLbAe2Jdfful6jY0aNcrGjh3rnlfybtCgQS6WDPocgLKjnKfiOshJusDXgVxdd0vi9sIAAOSsCRNU1dhs/Hh1VUjbZDnXlj6sMwAAMmjCBFu+446mW0yE6VxLzS8AAJD7VOxe9Wj0CAAAgMxo29bs1VctbBj2CAAAcp9aHQ84INtzAQAAEP6Yq3dvCxt6fgEAgNynwslXX533CAAAgMxQrHXTTRY2JL8AAEDpCMSuuYbkFwAAQCb99ZfZzTdb2JD8AgAAAAAAQGiR/AIAAAAAAEBoUfAeQFZs2LTZxk77x35fuNJWr9tk1atUsDYNato+HRtapQrk5QEAAEozb8MGWzlunK37Y6ZtXr3KylevYVVat7KavXpZuUqVsj17AMqYlJNfW2+9tc2ePTvy/+XKlbMaNWpYnTp1rF27drbjjjvawIEDrUePHoGf32uvvWzcuHH2ySefuH8DKFsWLF9rz307x579ZrYtWbXeKpQvZ+UUIJnZps2e1atR2Y7ZpaUdvXMLa1S7arZnF0Cu2GILs6OPznvMQcRHAJBnw4J/bOmLL9q/o0bZpiVLzCpU0EHRzPPMNm2yCvXq2RZHHml1Bw2ySo0aZnt2AcRSrDVwoNlLL1mYFLnnV8+ePa1t27bu32vWrLFFixbZjz/+aJ9++qmNGDHCevXqZY8//ri1bt06nfMLoBT7+vfFdtJT39uaDZtss7Jd/0t4RVNC7L6xM2zk53/YY8d1t13bbJmdmQWQW1q1Mnv2Wct1xEcAyrJV335nc08/3by1a802b857ctOmfO9RQmzRQw/Z4ieftOYPPmg1dg5uFACQxZjr0UdDl/wq53lKwafesvnEE0/Y8ccfn+81Terdd9+18847z2bMmGGNGjWyr7/+2lpp4f3PnDlzbPXq1daiRQurXr16+n5JCC1fvty1GC9btsxq166d7dkBip34Onbkt7bZ8yKJr0TKl9NfOXvmxJ1JgAEw04XUn3+aNWtmVrVqzp1riY9KDvERkLuJrzknnpiX9PITX4mUL+/+WowcSQIMyCVr19ryX36xOjvsEKpzbVoL66iLf58+fey7775zXfwXLFhgJ510Ur73KKjr2LEjgR1QxoY6qsdXsokv0fv0/pOe/t59HkAZN3WqWbt2eY+lDPERgLIw1FE9vpJOfMn/3vvnGWe4zwPIEVOnmu2wg4VNRqpK161b1+666y7377Fjx9r48eMjr6mOhYJAdf/37brrru65UaNGxZ3mfffd597Tr1+/Aq/pc/vuu6/Vq1fPqlSpYi1btrShQ4far7/+Grd1VtOaNWuWvfHGG7bPPvu4z8bO17///mvXXnut7bTTTq6FsVq1am6Ygmp2qAU31saNG+2xxx5zv9GfF7Xqnn766TZ37twUliAQLqrxFT3UMVl6/5r1m+z5b+dkatYAoMQQHxEfAWGlGl/5hjoma/Nm27xmjS0N2fAqALknY7dUO+igg1yAIx9++GHC955wwgnu8cknn4z7Hg0jEAVt0cMIjjvuOBs8eLB99tln1q1bNzv88MOtatWq7v36//feey/uNFV747DDDrMVK1bYgQce6OpwVFBBRjP76aefrHPnzjZ8+HD77bffbPfdd7e+ffvaVlttZW+//bbdcsst+aalaey333528sknu2B2++23t0MPPdQFeA899JCbF9X8AMriXR1V3D7VxJdPn9PnNR0AKO2Ij4iPgDDe1VHF7VNOfPk2b7Z/R73gpgMAOVfwvjBqJdxhhx3so48+silTpiR875FHHmnnn3++CwLnzZtnTZs2zff6pEmTbMKECa5GhoJG38MPP2xPP/201a9f3322a9eukaDvmmuucX8K/NTC2aBBgwLf++CDD7qWTQVh0VatWmWHHHKIm5chQ4bY/fffbzVr1oy8rnGv33//fb7PnHbaaa5V9D//+Y+NHDnSGjb8/zuXqJVXv2/QoEH2yy+/RALIbFq9fmO2ZwFlxEe/LHBF7Itj8ar19t7Pf9u+23BHIGRe9coZOzUCxEc5Hh9tXr0627MAlDorPvkk766OxbBp8RJb/uGHVou73QIpK0/JhKRkNMJX0CWLFy9O+D4VUDviiCPsmWeeccHaZZddFtiqeeyxx1rFiv8/y7fffrt7vOqqqyKBnR9YqkVy9OjRLjB89NFHbdiwYQW+V62isYGdqGu+uuFrmrojU2wwpi7+vXv3jvy/ArYXXnjBmjRpYs8//7zVqlUr3/tV4FbB55gxY9xwAAWAQdatW+f+ogu6Zkqnq97P2LSBTDj7BXoGoGTMuvngbM8CQo74KHfjo+k77JixaQNIbP4FF2Z7FoBSaZtpv2R7Fsr2sEfZ/L+urwq2CuN37X/qqafyPb9hwwZ77rnnCnTp//PPP+3333+PBGmx9J3+ND/55JPA7+zfv3/g8/5QgBNPPDGpVkgFbWpNVatrbGDnU50L+eqrr+JO56abbnKBo//XvHnzQr8bAIAyQYVXdYPqEBRgJT76f8RHAADkmB12UHduC5uM9vxatGiRe/RrWySi4EfFUqdPn+4CoN122809r/oRCxcutJ133tm22WabyPvV5V623HLLuLfebNOmTb73BhV2DaJblYvuupSMP/74wz2qO7/+EtFviUctuhdccEG+ls1MBXhTrz0gI9MFYj362R92z9jfbFNRi36ZWYVy5ezc3m3tpD1ap3XeACAbiI9yNz7qMOH/b0IAIDmLn3jCFj3woNmmTUWfSIUKVv/MM2zL449P56wBQOaTX2rp8wuYqjBqYdQSefzxx7su+irs6gd3fpd+v5UynXR3onS24GoYQJcuXRK+V0FqPCr+qr+SQE0blJSOjWsXK/ElmzzPOm5Vm+0WKMumTzfTRZGKv3foYKUV8VFux0fUTQFSV1XH5OIkvmTTJjcd9kEgR2KuY46xsMnYlaS6uutW2LL//vsn9Rl1z7/66qvtxRdftLvvvtu17KkGhIIwFX2N5hd9Vb0MvS+oddNvcYwtEFuYFi1auDoV06ZNy1e7Ih6/9bFnz57uluMA/t8+HRtavRqVi1X0fssalW3vjhS7B8q0VavMvvkm77EUIz4CEDY1dUfYevWKVfS+wpb1rOaee6Z1vgAU0apVZj/8YGGTkZpfutuP7t4jur11dLHVwoKqfffd1wVrr732mj377LO2ceNGd3tu1XiI1qxZs0i3/aBbgKtl1X9+7733Tmn+dVtvUTHXTUm0Yvh3WHrzzTdt7dq1KX0XEHaVKpS3Y3ZpaeULL20TSJ/T5zUdACjNiI8AhFG5SpVsCyXiyxcxVitf3rY4crCbDgBkSlqvJhVQqSWyR48eNmPGDGvcuLG7k1Aq/KKt6s5fWJf+//73v+7xuuuus59++inffFx//fU2ceJEq1u3rp188skpzcNJJ53kgkcNS9BndWvvaAo+dYtyX7du3dzdmHQHJAWis2bNKjBNTUOFaRcsWJDSvABhcPTOLaxapQopJ8D0/mqVK9hRO7fI1KwBQMYRHxEfAWFXd9AgK1e1auoJsPLlrXy1alZ34MBMzRoAFG/Yo253/emnn7p/6/bTKt46YcIEW/K/7q4q0KqWwZYtW6Y03cMOO8y22GIL+/jjjyNFV/fZZ5/A95566qmu+KtuAb7TTjtZr169rGHDhm4+VBhWwwF0a+0GDRqkNA81a9Z0rZR9+vRxAaZuCa4u+3peAZyCPgWw0V3+9b6lS5e64LZDhw6utkWrVq1coKlgT8Hn+vXr3XCBRo0apTQ/QGnXqHZVe+y47nbsyG91+WXJlABT4qt8uXI28rju7vMAUBoQHxEfAWVRpUYNrfmDD9qcE0/Me+J/Nf8SUqKsfHlr9uAD7vMAkJPJry+//NL9SY0aNVy3exVuVZA1aNAg6969e5GmW7VqVRs8eLA98MADkToX8W4Frueffvpp163+kUcesfHjx7sWxK222soVh7300ktdoFUUaq2cPHmyq63xxhtvuEBWhVvVWnvooYcWaG3VLbw/+OADV49DwxE0L2pZVa0Nfeboo492n/OHIgBlza5ttrRnTtzZTnr6e1uzflPCBJjf40uJr11ab1mSswkgV+kOhM88k/eYw4iPiI+AsqrGzj2sxciR9ucZZ9jmNWsSJ8D+1+NLia8aPXqU5GwCKIxirUceMTvlFAuTcp6a3pCTNHxAQbNqhMS7XTlQ2ixYvtae/3aOPfPNbFcEv0L5cqbLNx2IdFdIFbdXjS8NdaTHF4BM41xb+rDOgNy2YcE/tvSll+zfF17IK4JfoYKy8hp77e7qqOL2qvGloY70+AJy0/IQnmtJfuWwMG5wgG/Dps32ybR/7PeFq2zVuo1Wo0pFa9OghrurI8XtARSwcKHZSy+ZqS5MisP1EuFcW/qwzoDSwduwwVZ+9pmt++MP27xqlZWvUcOqtG7t7upIcXsghy1caMufesrqXHRRqM61RR72CADFoQTX/ttule3ZAFBazJ1rdtZZZrvumtbkFwAgM5TgqrXvvu4PQCmLuS66yMKG7hUAAAAAAAAILZJfAAAAAAAACC2SXwAAAAAAAAgtkl8AACD31apltv/+eY8AAADIjFq1zPbZx8KGgvcAACD3tWtn9v772Z4LAACA8Mdco0eb1aljYULPLwAAkPs2bTJbvjzvEQAAAJmNuUKG5BcAAMh9P/2U1wKpRwAAAGTGTz+ZNW9uYUPyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFF8gsAAOS+zp3N/vkn7xEAAACZ0bmz2e+/W9hUzPYMAAAAFKpSJbMGDbI9FwAAAOGPuerXt7Ch5xcAAMh9aoE89NBQtkQCAADkjN9/Nxs0yMKG5BcAAMh9y5aZvfVW3iMAAAAyY9kys/fes7Ah+QUAAAAAAIDQIvkFAAAAAACA0CL5BQAAAAAAgNAi+QUAAHJf06ZmI0bkPQIAACAzmjY1u+EGCxuSXwAAIPc1amR2wQV5jwAAAMiMRo3MzjrLwobkFwAAyH3//mv28st5jwAAAMiMf/81Gz3awobkFwAAyH0zZ5oNHJj3CAAAgMyYOdPs+OMtbEh+AQAAAAAAILRIfgEAAAAAACC0SH4BAAAAAAAgtEh+AQCA3Fetmlm3bnmPAAAAyIxq1cy2397CpmK2ZwAAAKBQ22xjNmFCtucCAAAg/DHX55+b1aljYULPLwAAAAAAAIQWyS8AAJD7fvzRrEqVvEcAAABkxo8/mtWvb2FD8gsAAOQ+zzNbvz7vEQAAAJnheWYbNljYkPwCAAAAAABAaJH8AgAAAAAAQGiR/AIAAAAAAEBokfwCAACl47bbP/+c9wgAAIDM2GYbs2++sbCpmO0ZAAAAKFS1ambbbpvtuQAAAAh/zLVN+Bob6fkFAABy3+zZZiedlPcIAACAzJg92+yssyxsSH4BAIDct3ix2ciReY8AAADIjMWLzZ55xsKG5BcAAAAAAABCi+QXAAAAAAAAQovkFwAAAAAAAEKL5BcAAMh9jRqZXXpp3iMAAAAyo1Ejs/PPt7Ah+QUAAHJf06ZmN92U9wgAAIDMaNrU7OqrLWxIfgEAgNy3YoXZp5/mPQIAACAzVqww+/xzCxuSXwAAIPfNmGG29955jwAAAMiMGTPM/vMfCxuSXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAIPdVqpR39yE9AgAAIDMqVTJr3NjCpmK2ZwAAAKBQnTub/flntucCAAAg/DHXtGlmdepYmNDzCwAAAAAAAKFF8gsAAOS+yZPNmjXLewQAAEBmTJ5s1rGjhQ3JLwAAkPs2bDCbNy/vEQAAAJmxYYPZX39Z2JD8AgAAAAAAQGiR/AIAAAAAAEBokfwCAAAAAABAaJH8AgAAua9dO7NPPsl7BAAAQGa0a2f29tsWNhWzPQMAAACFqlXLbK+9sj0XAAAA4Y+59tjDwoaeXwAAIPfpTo+XXZb3CAAAgMyYN8/s6qstbEh+AQCA3LdggdnNN+c9AgAAIDMWLDC7804LG5JfAAAAAAAACC2SXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAg9225pdmJJ+Y9AgAAIDO23NLs2GMtbEh+AQCA3Neypdljj+U9AgAAIDNatjS77z4LG5JfAAAg961ZYzZlSt4jAAAAMmPNGrNffrGwIfkFAAByn4Kw7bYLZTAGAACQM375xWyXXSxsSH4BAAAAAAAgtEh+AQAAAAAAILRIfgEAAAAAACC0SH4BAIDcV66cWeXKeY8AAADIjHLlzCpVsrCpmO0ZAAAAKFS3bmbr1mV7LgAAAMIfcy1aZFanjoUJPb8AAAAAAAAQWiS/AABA6bjt9g475D0CAAAgM375xWyPPSxsSH4BAIDct2aN2Y8/5j0CAAAgM9asMZs0ycKG5BcAAAAAAABCi+QXAAAAAAAAQovkFwAAAAAAAEKL5BcAAMh9rVqZvfRS3iMAAAAyo1UrsyeftLAh+QUAAHLfFluYDRiQ9wgAAIDM2GILs379LGxIfgEAgNy3YIHZHXfkPQIAACAzFiwwu+8+CxuSXwAAIPfNm2d24YV5jwAAAMiMefPMLr/cwobkFwAAAAAAAEKL5BcAAAAAAABCi+QXAAAAAAAAQovkFwAAyH116pgdckjeIwAAADKjTh2zAw+0sKmY7RkAAAAoVJs2Zm++me25AAAACH/M9eKLoWtwpOcXAADIfRs2mC1cmPcIAACAzNiwwWzRIgsbkl8AACD3TZ5s1rBh3iMAAAAyY/LkvN5fIUPyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFVMdszAAAAUKguXcyWLTOrUSPbcwIAABDumGvuXLPmzS1MSH4BAIDcV6GCWe3a2Z4LAACAcKsQzpiLYY8AACD3zZhhdsABeY8AAADIjBkzzPr1s7Ah+QUAAHLfihVmH3yQ9wgAAIDMWLHCbOxYCxuSXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAIPfpdtv33Re6224DAADklObNzW67zcKG5BcAAMh9DRqYnXlm3iMAAAAyo0EDs1NOsbAh+QUAAHLfkiVmzz6b9wgAAIDMWLLE7MUXLWxIfgEAgNw3a5bZscfmPQIAACAzZs2i5xcAAAAAAABQmpD8AgAAAAAAQGiR/AIAAAAAAEBokfwCAAC5r0YNs112yXsEAABAZtSoYbbTThY2FbM9AwAAAIXq0MHs66+zPRcAAADhj7k+/tisTh0LE3p+AQAAAAAAILRIfgEAgNw3YYJZuXJ5jwAAAMiMCRNC1+tLSH4BAAAAAAAgtEh+AQAAAAAAILRIfgEAAAAAACC0SH4BAAAAAAAgtEh+AQCA3Nepk9mMGXmPAAAAyIxOnUJ5g6GK2Z4BAACAQlWtata2bbbnAgAAIPwxV5s2Fjb0/AIAALlv5kyzY47JewQAAEBmzJxpdvLJFjYkvwAAQO7791+z557LewQAAEBm/Puv2UsvWdiQ/AIAAAAAAEBokfwCAAAAAABAaFHwPod5nucely9fnu1ZAQAgu1au/P/HNJ4X/XOsf85F7iM+AgAgg1auNP8MG6b4iORXDluxYoV7bN68ebZnBQCA3NCrV8bOuXXq1MnItJFexEcAAJSMFSGKj8p5YUrlhczmzZtt/vz5VqtWLStXrly2Z6dMUquyguu5c+da7dq1sz07SCPWbXixbsMpU+tVYZACuyZNmlj58lSDKA2Ij0oWx9Syi3VftrH+y67l/1v3U6dOtQ4dOoQmPqLnVw7TRtasWbNszwbM3AGfg344sW7Di3UbTplYr2Fp0SwriI+yg2Nq2cW6L9tY/2VX06ZNQ5P4kvD8EgAAAAAAACAGyS8AAAAAAACEFskvIIEqVarY8OHD3SPChXUbXqzbcGK9AtnBvld2se7LNtZ/2VUlpOuegvcAAAAAAAAILXp+AQAAAAAAILRIfgEAAAAAACC0SH4BUVasWGHDhg2zDh06WLVq1ax+/fp28MEH29ixY4s93RtuuMF22mknq1Onjpt2y5YtrV+/fvbRRx+lbf5Rsus12sCBA61cuXLu79lnn03bdFFy63b16tX29ttv21lnnWVdunSxWrVqWeXKla158+Z25JFH2pdffpmR31BWvfzyy7bXXnvZFltsYTVq1HDL/NZbb7UNGzYUaXrjx4+3AQMGWKNGjaxq1arWqlUrO/vss+2ff/5J+7wDYZXu8+V9991ngwYNsk6dOrlpVapUyerVq2d77rmn3X///UXe35EZnFPLrnTv+xMmTLDbb7/dBg8ebO3bt7fy5csTI2cRMdf/qOYXAM9bsGCB1759e9XA8xo3buwNGDDA23PPPb1y5cq5v3vuuadI0508ebLXtGlTN91mzZp5hx12mNe/f3+vR48eXqVKlbwLL7ww7b8FmV+v0UaNGuWmr+np8ZlnnknLvKNk1+2jjz7qpqW/li1ben379vWOOOIIr23btpH1e/3112fs95Ql5557rlumFStW9Pbff3/v8MMP9+rWreue23333b3Vq1enNL2XX37ZTUuf7969uzdw4ECvdevW7v8bNWrkzZgxI2O/BQiLTJwvFf9o3+zSpYvXp08f78gjj/T22GOPyP668847eytWrMjI70FqOKeWXZnY97W+/fUf/UeMXPKIuf4fyS8g5iC97777eqtWrYo8/84773gVKlTwypcv7/30008pTfPvv//2GjZs6A4QDz74oLd58+Z8r//777/elClT0vYbUDLrNXYdb7nlll63bt28nj17cmIvxev2ySef9IYOHepNmDAh3/Pab0eMGBEJ3D799NO0/o6yZvTo0W451qxZ0xs/fnzk+YULF3qdO3d2r6XSKDBv3jyvevXq7nMPP/xw5PmNGzd6xxxzTCQ4iz3+Asj8+fLzzz8PTG7Nnj3b69Chg/u+Sy+9NC3zj+LhnFp2ZWLfv+mmm7xhw4Z5r7zyivf77797vXr1IkbOAmKu/Eh+AZ7nElDaWXWAnzVrVoHXTzzxRPe6WixTceyxx7rP3XjjjWmcW2R7vUY79NBDXQ++iRMncmIP2bqNpaBQ09S0UXQKirQcg1r8daGs16pUqeItXbo0qelddNFF7jO9e/cu8JouuuvUqeNef++999Iy/0AYZeOY+vTTT7tp7rDDDmmbJoqGc2rZVVLrnhg5O4i58qPmF2Bmo0ePdo89e/Z0tbhiHXXUUe7xrbfeSnpstMY8jxo1yo2bV70DhGO9Rnv66aftzTfftMsuu8yNnUd41m2Qbt26uce5c+emZXpl0bx58+z777/Pt46i7b777q4ezLp162zMmDEpbQtB06tZs6Ydeuih7t+vvfZaMeceCK9sHFMrVqzoHqtUqZKW6aHoOKeWXdlY9ygZxFwFkfwCzOzHH390jypIH8R/ftWqVTZjxoykpvnJJ5+4k4RO7iry+dVXX9nll19up556ql155ZU2bty4NP4ClNR6jT6hnHvuubbddtu59YrwrNt4/Ok0btw4LdMry+tNBa9VHDXRuvPfW1iB3t9++y3f54ozPaCsKuljqhoIb7nlFvdv/2IJ2cM5tezKxrpHySDmKiivyQUo42bOnOkeW7RoEfh67dq13d/y5cvde3XXosJMmjTJPW611VZ27LHHFri7yfXXX2+9e/e2l156yd15A6VjvfpOOukkdxJ44okn3N2LEJ51G2Ty5Mn2zjvvuH8fccQRxZpWWVbYehO1Qka/N5FZs2ZF/h1vmqlMDyirMn1Mff755+2DDz6wjRs32l9//eXu9KfeBscdd5xdcMEFafkNKDrOqWVXSa97lBxiroLo+QX8L5MtuvVrPOrKKTr4J2Px4sWRbsIvvPCCXX311e5AsGTJEtcVVC1dH330kbvdM0rPepVHH33U3nvvPbvooovitnygdK7bICtXrnTdu3XRdsABB9ghhxxSrOmVZeleb/70Ek0zXdsBEGaZPqZ+99139tRTT9lzzz1nY8eOtfXr19v5559vd999Nw1IOYBzatlVkuseJYuYqyB6fqHUu/jii13dpVQ99thjbqxzpuiGEqKhj6oJNXz48Mhr/fr1syZNmtiuu+7qWkI///xz22OPPTI2L6VRrq7X2bNn24UXXmjbbLONS2giPOs2iPbfAQMG2M8//2ytW7e2Z555pkS/HwDCcEy966673J96e6n3gHrDjxgxwjUGqtYMvUnCvf59nFPL7roHcgHJL5R68+fPt+nTpxep5cmnmlz+ePbC3q+uv8nwpymq8xVr5513th122MHGjx/veoCR/Mr99aqE5tChQ930Hn/8cYr0hmjdBlGrtHpmqpefisCqt0KDBg2KNC1kZr1FH2c1zTp16hRrekBpVFqOqaLzZocOHey6666z7bff3gYOHOiGPvpFmRHe9c85teyue2QHMVdBDHtEqafWQyUlUv078MADI9PYeuut3eOcOXMCv0NdN/3um/57C6MWLf9uRv7453jvUf0L5P56XbZsmQvWqlevbpdeeqnttdde+f4mTpzo3nfDDTe4/z/vvPPStDTCJRfXbaxNmzbZ0Ucf7XolaP/VDSyC7oKE1PjrItHdvfzXkllv0esk3raQyvSA0qg0HFODqNaTLqZ++OEH7vgX8vXPObXsrntkDzFXQSS/ADPXA0sUgAXxn9f45vbt2yc1zR133DHS0hVv3POiRYvyjY9G7q9Xv1VDd+uM/VNyTKZNm+b+30+GofSsWz9IP+aYY9zNKPwgPd5dclC0W9urJmK8Yqj+uvPXcSJqWWzbtm2+zxVnekBZlcljajzly5e3atWqRe7+iOzhnFp2ZWPfR8kg5iqI5BdgZocddph71N2HgjLZukuRqChnpUqVkppmjx49InfC+PDDDwu8rsL3GvLovxe5v17r1q2bsCWtV69e7n2qYaH///TTT9P+m5C5fVY2b95sQ4YMsVGjRkWC9DZt2qRxzsu2Zs2aWffu3fOto2hffPGFazXU0Kg+ffokNU3VUIw3PSWqddMROfzww4s590B4ZeqYmojqPinpVaFChUhPeGQH59SyKxv7PkoGMVcAD4DTt29fVaj3evfu7a1evTry/JgxY7wKFSp45cuX93766acCnzv22GO9Dh06ePfee2+B10aOHOmmufXWW3vTpk2LPL9q1Sqvf//+7rUWLVp4a9asyeAvK9sysV7j6dWrl/uuZ555Jm3zj5Jbt5s2bfKGDBniptm8eXPvt99+K5HfUdaMHj3aLeOaNWt648ePjzy/aNEir3Pnzu61Cy+8MN9nXnvtNbfO9tlnnwLTmzdvnle9enX3uUceeSTy/MaNG9261vPdu3f3Nm/enOFfBpRu6T6mfv75596bb77pbdiwocBntO9vt9127vuOPPLIDP0ipIJzatlVErEyMXJ2EHPlR8F74H8eeeQRmzp1qis+r1YpFaBXi6SGr6kXj27HreKssdRKomKT/hDGaCqO/vXXX7u7qnTt2tV22WUXVxzw22+/tb///tvq1atnr7zyilWtWrWEfmXZk4n1inCu2/vuu8+efvpp929NTwWZg3Ts2NHVfEPRW5nPOeccu+eee9wxcd9993XDKT7++GNbunSp9ezZs8Cy15BirbO1a9cWmJ7unPvkk0/a4MGD7ZRTTrGRI0e6WhMqoP3HH39Yo0aNXAtluXLlSvBXAqVPuo+pv/32m51wwgmu17SG3zRu3Nj1DNDwm8mTJ7v3aH9/8MEHS+w3Ij7OqWVXJmLld955J9861/RFd0rXtuH75ptvMva7QMxVQEwyDCjTli1b5l166aVeu3btvCpVqnj16tXzDjzwQO+jjz4qtCVj+PDhcd/zwgsvuPfVqVPHq1y5stemTRvvrLPO8ubOnZuhX4KSWK/xPkOrVulct/p/PV/Ynz6P4nvxxRe9Pffc06tdu7ZXrVo11wvk5ptv9tatW1fgvU888YRb9i1btow7vR9++ME7/PDDvQYNGrjjrN575plnen///XeGfwkQHuk8ps6cOdO78sorvb333tv1/KlatarbN5s1a+Ydcsgh3nPPPed6ByF3cE4tu9IdK/vn7cL+UDKIufKU038KpsQAAAAAAACA0o+C9wAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwAAAAAAAAgtkl8AAAAAAAAILZJfAAAAAAAACC2SXwBKzIwZM6xfv37WuHFjK1++vNWtW7fQ166++morV66cffrpp0X6Tn1On9d0SpPXX3/dzfdXX32V7VnJadOnT7eKFSvaAw88kO1ZAQCE0F577eXOx8l68skn3fv1WFbcc889tu2221r16tXdb7/rrrsKfU3/1rItqfWSCzzPsx133NH233//bM9KzjvmmGOsZcuWtnbt2mzPCkKE5BcQYuPHj7cTTzzR2rVrZzVq1LBq1apZmzZt7Nhjj7UPP/ywROdl06ZNdthhh9mYMWPs4IMPtquuusouvfTSQl8rrY4//ngXlM2aNSvlz27YsMEuvvhiO+CAA2y33XbL99rmzZvtvvvusx122MEFkrVr17Y999zT3nzzzcBp+cnDeH+x86d1ceONN1rr1q2tTp06tt9++9nUqVMDp/3II4+4efj999+tqEHga6+9Zocffrg1a9bMqlSpYrVq1bIuXbrY+eefX+B7/UD377//jjzXoUMHGzx4sF1zzTW2YsWKIs0HAKBwOl/oGHzggQdme1bKtNWrV9vdd99te++9tzVo0MAqVapk9erVs913391uvvlmW7hwYYnOz6hRo+zcc89153A9Dh8+3HbZZZdCXyvN+4BivKJ4+umnbcKECXbttdcWeO2XX36xo48+2rbaaiu3vJT40TJbsmRJ4LQSxXZB86eYqnfv3i6207XATTfd5GK+WGvWrLG2bdvaKaecYkU1b948u+yyy1ysqsbsypUru8ZtxfhKCq9fv75AI/Vpp52Wbxq6FtB0ohOpQHFVLPYUAOQcJUj++9//2p133ul6xeyzzz526KGHugDpjz/+sHfeeceeffZZd/K98sorS2SeZs6c6U68J598skuaJPvaWWedZUceeaS1aNGiSN/bo0cPF1DUr1/fSotnnnnG9YR76KGHCiSLBg4caK+++qoLXJTYXLdunb3xxhvWt29fu/fee93yCnLcccfZ1ltvXeD56N53om3m8ssvtz59+lj79u3ddqIWymnTplnNmjUj7/vrr79cgk7JNc1LqhTMDRgwwMaOHevmQUk2JdwUEE2ZMsX15FJr8ccff1xoy7DmQ/Op92veAQAIo59++smd72fPnu2SI4rtGjVqZMuXL7dvvvnGJRyU1Jg/f75r9CwJb7/9duSxSZMmSb+m2EwNaEWlRJISgaUpNlfMtMceexRIAGrdKTGlxJPWr+KqiRMnurjmvffec6MAttxyywLT1DYQlOjq2rVrvv9X46Cmv3HjRhc7atkPGzbMJdkuuOCCfO9VglLL9bbbbivS73zhhRfcd+i3qJebenAp4aaGS8V8J5xwgotzFd8lohhUy0IJ3bPPPrvEtmeEnAcgdC677DJPu3fXrl293377rcDrq1ev9m699VbvkksuKbF5GjdunJun4cOHp/RaaXXccce53zRz5syUP7vTTjt5zZs39zZv3pzv+ZdfftlNs2fPnm4d+hYuXOi1bNnSq1KlSoHv0zLVZz755JOkvrt9+/Ze7969I///2Wefuc+/8MIL+d7Xr18/r1u3bt6GDRtS/n36zJ577umme8wxx3jLli0r8J758+d7J5xwgvf6669HnuvVq5f7zF9//VXg/dtvv71bBps2bUp5fgAAhdP5RcfgAw44wCtL/HNPsp544gn3fj2m09y5c72GDRt65cuX90aMGOFt3LixwHsmTJjg7bbbbt6///7rlZS999477vJJ9Fpp3gcU46Xq7bffdp999NFHC7y23XbbudfeeOONfM8rVtfzp556aoHP6Hltm8l4/vnn3fu/+OKLyHP77LOP16FDhwLbT8WKFb1XX33VK4p3333XbZ/16tXzPvjggwKvK6597bXXvD59+kSeU3wa7zfqvXrtscceK9L8ALHCczQC4MyYMcOrUKGCt+WWW3p///13wveuXbs23/8riXLuued6W2+9tVe5cmWvQYMG3oABA7zJkycHfn7dunUuAFMSpHr16l7NmjW93XffvcDJW0kJnbxi/5SYSfRaYcmbiRMnekcddZTXtGlTN79bbbWVC8rffPPNAifVoMTaggULvPPOO89r06aN+7yW2eGHHx74ezWf+luxYoV3zjnneI0bN3af6dy5s0tKJfN7kwlS9N16r+YrlhJFeu2dd94p8Npdd93lXrvqqquKlfyqVq1avqSokmz6/C233BJ5TkGRtrHx48d7RfH444+7aSoBVliyKnobTZT8uv76691rH330UZHmCQCQ3uTXrFmzvKFDh3pNmjTxKlWq5M7V+v/Zs2cHvv+nn37yDjroIBdL1K5d2/1b58SiNibps4phFMvofK3YRjHOokWLAt//+eefu/OS4hldvA8cONCbM2dO3OTX4sWL3QW7ElI6d6rhShfr8ZJfY8eO9Q488MBI/KDPKWZ6+OGHk/o9Q4YMcdO94oorCm1gij23Ki7aa6+93HKtWrWqazBS/BavAUvrYtCgQS6u0rpr0aKFd9ZZZ+Vbdv7vDPpL9JovXlyk2PKOO+5wy1PbQo0aNbxtttnGO//8870lS5YklZRUw5mSO3Xr1nUNg9tuu6132223FUgYRq+r999/39t1113dutT61/JO9vcmE2MdccQRXrly5fL9BlEjtabRvXv3Ap/RelRsqm1y5cqVRU5+KYbT+9esWRN57uKLL3a/1adlo3hejZtFoc+3bt06qVgsOrZLlPzStqDfrv0ESAeGPQIho7H0GsN/6qmnuq7wiai7s081InbddVdXv0nDzDTUUMMRX3nlFTdM8v3333f1JHwabqe6Hxqrr+7V6uKsWlV6b+wQvPPOO891337qqaesV69ekWFsetSQt3ivJaKhf0cddZQbCnjIIYe42k///POPffvttzZy5Ej3XCL+7/zzzz/dsD7VHNPnNV39VnXH3nnnnfN9Rr9P7/3333/tiCOOcN3CVc9CQxHVLd0vYKrfq/Wg4Qmq1+APLQwadhjL7wYeVBPDr3XVqlWrAq/5z6lLuepfxfrss8/cstHNBFQDTt3fo4cx+po3b24//vhj5P9Vm0L8YafLli1z61U1uVTLoSi0fuSKK65w85PsNpqItl1/+e27775Fmi8AQHr8+uuvLmZQbKHzsQqe//zzz/b444/bW2+9ZV988YUb1uTT+VLDwVatWuXqQOo89cMPP7hpqA5kqjR91c3UUPr+/fu78+/XX3/tamVpCJ6GmUWXQ9C546CDDnLnpEGDBrkhenquZ8+etsUWWxSYvs7/iiEmT57szj+KX+bOnes+G1TMXLGRloPiAcVIqn+kZaPfrSFghdVX8uMN1W5VWYtEVO4i2h133GEXXnihqwumuEnDx1QnVM99/vnnrvZmdOF4vaa4RstC86q4QKUpVG9U8ZFiCS0TxX4aIqd4R8Mw9W9fotcS0VA5lUH48ssv3TagIXKKA1QK4uGHH7YhQ4YEro9oGvqpoXJNmzZ125KG3Ol3XnTRRW7eX3755QKf0W/215FqrSpm0rBKxYralvzfpJhO25C2ScWNvsLiO8Wqn3zyiYtVY+c/UWyndaD4S3GZttnY+Gbp0qWuXMiiRYvc+tX22rlz5wLT0ToUTcePlxTfRZcUGTFihCuN4g9VTZV+nz6v5VdYHJZsbKdaYRo6qX1XxwaGPqLY0pJCA5Az1LJXlB4wGmKmz2nIZDT1MtLzbdu2zdeSOGzYMPf8lVdemW943vLly11rnVo1582bl1QPrESvBfVcUo82tQTqT120g4YGFDZtDQtQ76X33nsv3/PTp0/3atWq5Xp0BfXm6tu3r2uJ8mk5B7WEF7WlWq3U+px68MU68sgjC+35pZbkoOUX+6fW0KeeeqrAdPwu9gcffLBrZVWLuVrt1eNNTj75ZNeyFz3sMhVqZVYrsrrVR7dAJiNRzy8NnfR7kwEAstvzyx/uFtur6f7773fPq1dONPXs0PPPPfdcvucVY/jnrWTPp4pV1KNbn4k9x1900UXuefVAi36/zmvqlaPeXz7FNupdHttjKfrcqnNiNH1fdO8nn3qV6zn1WI8VrydatE8//dR9PtUeMOpVpPOtYgP1YovueeMv86effjrfvKh3mHrpqedeNJU/0PvVAyxaoh5YiV4L6rl04YUXuuePPfbYAr20li5dGolF4k1bQ+38bTS6p5TW5WmnneZee+WVVwr05tIyih4SqO/24+mvv/662MMep0yZ4j539NFHF3ht2rRphfb80usPPPBAvtfi9UJT70KNbIiNkdSLr1GjRi6207BDvVe9//ztRL3AHnroIa+orr766qR6JsZK1PNLNL96XT0ngeLibo9AyPgtSLp7XrLUMqoClSqmqd440VT4XK1wv/32m2uJ84t2Pvjgg64gp3oZRbcY6m59ukOLpqnWxExQLzG1AKnVslu3bgVeL+y3q+VLxUNVBF4tw9HUEq3C+2rNVSt1LBWEV0uUT61bKjj6/fffWzqoJ5oE9dpTq7SoRTP61s+LFy+O3A1HrYDR1Dqplna1xqlFVb351CvPvxtQ7F0iVfj0uuuuc0Xn1UNr++23d3cGVS+xcePG2WOPPeZaX3XzBL8luWrVqpGec4XRvKoHnVrc9bl00V0vNT1/+QEAsmPOnDmuF0inTp3c+TSa7ujWsWNH10tZPaVEPYPUu0bnK/VMinbJJZcU2tMnlmIV9djROTP2HK/4ROet559/PnLHOX23zpH/+c9/8vVw13lSdz+uUKFCge9QryDFArF37dP3Jer1op5bsYIKmacjthP9ThU51/na7/3j97y55ZZb3L/VOyv6d6l4vormK7aJphEB6vGtHmiZoPlULyb11FLvqtjlrueDeqxHU+800XSiewlpXSp20qPi3Vja7tRryqfvVowo6YjvEsV2ijt1wx99j3qfRVNsp7gpKL7TOlUsq15fWmf6t7Z5jUTQthx9J0fFSIrl1ANTcZwK3t9www2uJ5totMhOO+3keiCqd5z+rR6E2gYUcyejqNtoYfxlRnyHdGDYIwB3Jz8lU3Tb7KA77+h5nTQ1PFHDEqZPn+6G/mlYQNAQO/8225puJnz33XfuMWhoQTLUdVwWLFjg7rwTy59vPW633XaR5zVcIahbuk706pKdDgpyFHQpiRgUnClI1UWFurVr2KkSSa+//nokOIgdRtivX798/6+u+Rq2uM0227ikppKduluUT9+t52KToNo+FBRpyIGGTCowvv/++10yUL9f01QybcyYMZYtuqBREAgAyB7FCqKhgNGNY/45as8993TnV71PCRkN/ZPo5INPCQwNN9N5z6ckgN/gE80/n/tD94PKJyh5ogv7Dz74wMUyOpf636/4JpYu/jWPs2bNijynRIMakpTc22qrrQp8RtOJvZOdEkdqEFRJA53LlSDT+zJ9J+pEy0LD39Ro5K+v6PhICRAlEGMpFtB5Vn/pnndtE/5dCVNNeEbPv7YZNfrFSz4GxaYaWhfLT+LEJp2Kwk9gxd5hW7SP6A7XGnKpeEyNiWpc1nap7VTbqBpkY+O722+/vcD61JBF3eFdjZW6E7iGffoUzwbdYVHLSglgfd/KlSvt4IMPdvuckmianmI7JaxjS4GUZGwnxHdIB5JfQMgoENOJfd68ea62QDIUyEm8GmGqTRH9viVLlrhH9Q7SXzzqnZUJqjslqudQFP78q4UttpUt0fyr1TGIWsfUGy4dFJiptU5JLfWuiv2ed99917VeqjXXbyFVgks1QNR62LBhw6S+R4G3gisFVFqvahVMRK3bCgBVO0TUKqtE2Omnn+7+X9PQ7ax1MZFou1MLt36XAkHVjUu27kMy1LOtOLdNBwAUX6oxhf8Y7/wVOx2di4Ia3vzkV6rf78cUib4/NvmVyvzKgAEDXEOVzqEPPfSQazxS0kONi6q1pGRDIn6STbFdKhItC32/no+eph8faf4SUXyU7uRXcWM7f/7Vgyxo+0gUmwbFQH7ttOgeVEXl9/iL7rUf22NQdcnU8169IhWbKlk1evRol7BSrJZMfKcEmXpbKvmlHpDRya8gagRW/KgGT8Vu2ja1DNXQquSfEpFq1FSyOajHXDq20WRiOyG+Qzow7BEIGb/lNKh1Jx7/pK+TYKKuzP77/EcVff/fXWMD/5544gnLBL/lrKgnWH/+Nfwv0fz7Xd5LUoMGDfIFoLGULFLhWCWZlDzSUEMNQ/SXhVq0k+UHriqkm8ikSZPstttucwkvtcApQNW2Eh2s+8NPC+vtp2CyR48eLrmngrLpouSj5stffgCA7ChqTBFv6HzsdNSDOeicXdTv9xu2kv3+VOfXp+LxSkqo57wask466SR30yD14i6sd1H37t3dMEvdBMBPaCUj0bLQMtPz0Ykf/99KtiSKj2KHROZCbOfPvxrZEs27eu3lWmwn6lmlnlbaPpQk07pWLzCti1TiOz+2S6YB+uyzz3ZJLg0vFsWW+nz00EXFesmM5CjK9Ucy/GVGfId0IPkFhIy6J2vomnoF+cMP41HyRNSdWV3fVW8gKBGi4Ez8ZIeGzCnA0IlZSYySpuSJqDt4Ufhdt9M1VDGIX6si1RZD/y49CkBS8dxzz0WGViRDQZF67Wl4QKLWW82/AnS1SsZO299+ov8dO8QliO4MKqqlEn3BEiT6OxLRnaCUAAu6yxEAoOT4sYIaOGKP8fp/v+HDf59/N0fVLIqlmMQflpgsvzHGj11iz32KXdQTx++l7H+/et7EUj0yvzaZT/GPSiCoFqqfSIsWNJ1oKmughJfiNMVsSkBpmGEi6vWic7B6wainWCLq+eT3Rk+0LPSdSrJEN2SVRHwUj9aHlq1iUSWAikLzr57liglyKbZTrS31yko1tvPr4WmIbbLxjb8tFXYHSt11VXc4f/TRR/ONNIiNu/T/ycR26sWo2mXaj6OHKRcnthN/mRHfIR1IfgEh07ZtW7v44ovd2HgVvgxq4VKwo673/hABtSYOHjzYfUa1nKJpzL9uba3p+q066r2j4W46Kau7dFACTMXikymAXhTqkaW6HQoAo2tV+AprNVTyTAGSunC/+OKLBV5X0KjW2XTUKIgNmgujGikSLxAOavF95ZVXXM0GtQxHd3FX7Qzdbj6Wgmd1i9frup157G3Ro6m3lwqj6gYHPrWSa9hIdH0v/99KjBbm2GOPdbVOFIzrNuaaj1i6GNA8avtLhr+8/OUHAMiOFi1auAthNbDE1l5SwkfnFNUl8guwqxeR4gudz2PPyep1nKi3TBBNS8P61bvqo48+yvfa9ddf75Ijinn8m9eoyL2SWep1o0RDdKJu2LBhgYkOncdUMF8F9KOpUS6o54sSfkHT8eOkZG4AowLl6v2ix3vuuSew3IJ6aqu+lx8rqL6YzvGK+ebPnx95n+bd7+2jBJxP52Ql5y6//PLAshZKRvp1wdJN86nC6+rFrULssctLz6smVSLnnHOOexw6dGikzlY0JSu1/RWVapEpEZRqbKdebbqBkBKvQetNvys2Uazfq+1MyyE2NldvsKDYW4kn3chAySwNtY1H28cZZ5zh6rVG1/JSDKfX/BtcKT5TMjeZ2E6JQQ2XVZJPsaWGb8ZLuvXv39+SpfhOMWe7du2S/gwQDzW/gBBScKcEl4qRqyVNQaZqB+hkqGSYgkEFBXqfTydLJXz0nE6eOhmqxsXLL7/sWhw1hDG62KbqKUyYMMEFYKpNoAK2qkegxJNOymqpVcthsjWoUqFp6o5EagVVIksFQvU7lbzTSVKtXaqtkYgSXwrONQ3VMtAdjNQSrLtUab7Vay5ebYZkaJmrGKmKxGt4qHpYKcBXIFNYLS4FnrrBwEUXXVTgda0XXTAoEFGwrOL/SiKptU3rKvruSFrH6tWnpJjer3oMSipp/euuOWpF04VFPFr/CuwVdEXfJUoUmF566aWuMKq6x+sCR3cXSiY4UYCr9aPATHcR0h0ndfMCXXwoIJ86dar7TQrsVEcsGVpemq7mAQCQOTrHRydMoumco3ODGkyUVFIjhi521XNFyRQd75XAiW5Q8csQKI44+uijXW8UNbgpxlCiRc8reRRb8DsevU81i9RjWXes1rlG51+d23VuUWJMtTOj36+knN6rGkeDBg1yN/TRxftff/3lkhZKKkVTI6MK2KvXjH6X5lEJkZdeesmdF2PriSopo+STloliFCVQlGjTOVxF8KPvMhmPzrVKrmkonM7BivEUM6hulxIWmpZ6Tan3lN+TR79V8Z3uDKjfoaSE4hGtE/Wo0VDM6POs1o3iIy0z9YhTDzWtU/XUUUygOHG33XZLumEqVaovqnX+zDPPuEc14qrcg+7Gqe/UMktUH03ze+WVV7raWdqG9P9a94qH1FNPiRzFuckkc4Ko4VUxlbZHxXOKebT96N+FDQVVfVaVrdDv0jKMpphIiVbFjtr2lBTVvqJYVL8l+sZEosZfbWPabhSfaX1rO9T2oW1LSSit+3i0j2q+lUiNpmSp6n+pIVUJYu0vGpJ73nnnJbV8tLy17jRiQNumhmqqEL/iWsWfmp5upKD9LBl6r65b/PqyQLF5AELr+++/94YOHeq1bdvWq1atmlelShVv66239o466ijvww8/LPD+hQsXeuecc47XsmVLr1KlSl79+vW9/v37e5MnTw6c/saNG72HH37Y69mzp1e7dm03/RYtWngHHnig9+CDD3orV66MvPeTTz5Rk5Y3fPjwAtNJ9Jqe02t6T6wff/zRGzhwoNeoUSM3v40bN/YOOugg7+23305q2kuWLPGuuOIKb7vttnPLp2bNml67du3c8nnttdfyvVfLRH9BevXq5b4j1q233uqmp3nT63pfMk4//XSvQoUK3vz58wOXR+fOnb1atWp5VatW9bbZZhv3G5YtW1bgvXruzDPP9Lp37+41aNDAq1ixovtcjx493LytXr064Xzsv//+3i677OJt2rQpcN1fdNFFbhvRsjv88MPd9pOKzZs3e6+88op32GGHeU2aNPEqV67sVa9e3a0PbYdTp04NXM5//fVXvudXrVrl1p2mAwDIjJkzZ7pjcKK/6PPcrFmzvBNOOMGdm3X+0aP+X88H0Tn9gAMOcMdznat0Plf88Z///MdN+99//01pfidNmuRiGJ2ndB7WOfzcc8+Ne6767LPPvD333NOd0+rVq+cNGDDAmz17dtxz/OLFi71TTjnFnV91Pt5xxx1d7PDEE0+49+vRN2rUKBevtGnTxp3n6tSp43Xp0sW75ZZbvBUrVqT0u3TOu+uuu9x86bdp2datW9fbddddvRtuuMFbtGhRgc+88cYb7v1arorVFEeMGDHC27BhQ+B3TJs2zTvxxBPdMtO5eYsttnCf0bn5u+++y/feeMunsNfixUVr1671br/9dq9r166R2KxTp07ehRdemG8bSDRtxbiHHHKIWzda91tttZVbPtddd503Z86cyPuC1lVh8eP06dO9Pn36uGVerly5uDFqrHnz5rl1pRgv1sSJE912rn3Ej7/1/2PHjg2clrazvn37eq1atfJq1KjhPtO8eXNv8ODB3rfffptwPr744gs332PGjAl8XetXcaLWu6b/7LPPeqn6888/vUsuucTr1q2buz7Q71asrusDLev169cXWM6nnnpqgelcffXV7jUtHyAdyuk/xU+hAQDSRa2x6qmnYakaeoDEHnvsscjdjdT6DgAIBw35Ug8WDdePV0geKC3UQ0w9tlQ2RL2hkLh2nXrWaVRAvCGUQKqo+QUAOUZDONVlXEMaguphIX9wpML5GhJA4gsASu+xXKULYml4ohIFGuoHlHYacqlErob5IjGVxdC+rxIiQLpQ8wsAcpBqqqmOh2pscIeb+FSjbciQIYXWUgMA5C4V/G7atKntt99+1r59e1fzUTU8VcNKxa79G/QApZnqgimpQy/Gwql2mWrqqSYvkC4MewQAAACQNbrZiYpqa3iTCsPrhjNKeqnguQqYKzEGAEBxkPwCAAAAAABAaFHzCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAoUXyCwAAAAAAAKFF8gsAAAAAAAChRfILAAAAAAAAFlb/B97biLx1wpBEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ---------- 4) coef plots (1x2) ----------\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "ci = 1.96\n",
    "\n",
    "# left: Panel FE\n",
    "axes[0].errorbar(b1, 0, xerr=ci*se1, fmt='o', markersize=12, color='tab:blue')\n",
    "axes[0].axvline(0, color='red', linestyle='--', linewidth=1)\n",
    "axes[0].set_yticks([0])\n",
    "axes[0].set_yticklabels(['Divorce'], fontsize=16)\n",
    "axes[0].tick_params(axis='x', labelsize=16)\n",
    "axes[0].set_title(\"(a) Panel FE (Two-way): DV = Political orientation \\n(higher = conservative)\", fontsize=16, pad=10)\n",
    "axes[0].set_xlabel(\"Coefficient (95% CI)\", fontsize=14)\n",
    "if np.isfinite(b1):\n",
    "    axes[0].text(b1, 0.005, f\"{b1:.3f}\", ha='center', va='bottom', fontsize=16)\n",
    "\n",
    "# right: Cross-sectional Logit\n",
    "axes[1].errorbar(b2, 0, xerr=ci*se2, fmt='o', markersize=12, color='tab:red')\n",
    "axes[1].axvline(0, color='red', linestyle='--', linewidth=1)\n",
    "axes[1].set_yticks([0])\n",
    "axes[1].set_yticklabels(['Divorce'], fontsize=16)\n",
    "axes[1].tick_params(axis='x', labelsize=16)\n",
    "axes[1].set_title(\"(b) Cross-sectional Logit: DV = Conservative Vote \\n(1 = conservative)\", fontsize=16, pad=10)\n",
    "axes[1].set_xlabel(\"Log-odds Coefficient (95% CI)\", fontsize=14)\n",
    "if np.isfinite(b2):\n",
    "    axes[1].text(b2, 0.005, f\"{b2:.3f}\", ha='center', va='bottom', fontsize=16)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"figure_1.png\", dpi=600, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a967be-d5e2-4aa3-93a0-2bcecdafe1ec",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3062c375-e0e1-44f6-abf3-732ea6c39cc6",
   "metadata": {},
   "source": [
    "# 2. Analysis 1 and 2 full results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c494671a-c956-4be6-a9d8-2c8d65f1f419",
   "metadata": {},
   "outputs": [],
   "source": [
    "panel_df = pd.read_excel(r\"\\df_analysis_1.xlsx\")   # 복지패널\n",
    "cross_df = pd.read_excel(r\"\\df_analysis_2.xlsx\")      # KGSS 횡단면"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e95f80c2-f81f-4493-b492-27ee4232b116",
   "metadata": {},
   "outputs": [],
   "source": [
    "def safe_log(s):\n",
    "    s = pd.to_numeric(s, errors=\"coerce\")\n",
    "    out = np.where(s > 0, np.log(s), np.nan)\n",
    "    return pd.Series(out, index=s.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "fabf3b34-a273-4fad-be32-233e046c7c24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== Two-Way Fixed Effects Model Results ===\n",
      "                          PanelOLS Estimation Summary                           \n",
      "================================================================================\n",
      "Dep. Variable:               ideology   R-squared:                        0.0044\n",
      "Estimator:                   PanelOLS   R-squared (Between):              0.2524\n",
      "No. Observations:                9436   R-squared (Within):               0.0034\n",
      "Date:                Wed, Sep 17 2025   R-squared (Overall):              0.2510\n",
      "Time:                        16:10:48   Log-likelihood                   -8657.9\n",
      "Cov. Estimator:                Robust                                           \n",
      "                                        F-statistic:                      3.3303\n",
      "Entities:                        4145   P-value                           0.0015\n",
      "Avg Obs:                       2.2765   Distribution:                  F(7,5280)\n",
      "Min Obs:                       1.0000                                           \n",
      "Max Obs:                       5.0000   F-statistic (robust):             2.2708\n",
      "                                        P-value                           0.0262\n",
      "Time periods:                       5   Distribution:                  F(7,5280)\n",
      "Avg Obs:                       1887.2                                           \n",
      "Min Obs:                       1263.0                                           \n",
      "Max Obs:                       2642.0                                           \n",
      "                                                                                \n",
      "                               Parameter Estimates                                \n",
      "==================================================================================\n",
      "                Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "----------------------------------------------------------------------------------\n",
      "divorce           -0.4505     0.1412    -3.1914     0.0014     -0.7273     -0.1738\n",
      "ln_income          0.0094     0.0339     0.2772     0.7816     -0.0571      0.0759\n",
      "homeownership     -0.0839     0.0539    -1.5572     0.1195     -0.1896      0.0217\n",
      "ln_house_price    -0.0047     0.0224    -0.2112     0.8328     -0.0486      0.0392\n",
      "education          0.1228     0.1539     0.7975     0.4252     -0.1790      0.4245\n",
      "age                0.0491     0.0396     1.2419     0.2143     -0.0284      0.1267\n",
      "religion           0.0008     0.0413     0.0190     0.9848     -0.0801      0.0817\n",
      "==================================================================================\n",
      "\n",
      "F-test for Poolability: 1.7140\n",
      "P-value: 0.0000\n",
      "Distribution: F(4148,5280)\n",
      "\n",
      "Included effects: Entity, Time\n",
      "Optimization terminated successfully.\n",
      "         Current function value: 0.635028\n",
      "         Iterations 5\n",
      "\n",
      "=== Logistic Regression Model Results ===\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:           conservative   No. Observations:                19189\n",
      "Model:                          Logit   Df Residuals:                    19168\n",
      "Method:                           MLE   Df Model:                           20\n",
      "Date:                Wed, 17 Sep 2025   Pseudo R-squ.:                 0.01942\n",
      "Time:                        16:10:49   Log-Likelihood:                -12186.\n",
      "converged:                       True   LL-Null:                       -12427.\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.272e-89\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         -1.1313      0.104    -10.881      0.000      -1.335      -0.928\n",
      "divorce       -0.1857      0.086     -2.153      0.031      -0.355      -0.017\n",
      "income         0.0055      0.003      1.584      0.113      -0.001       0.012\n",
      "education     -0.0047      0.015     -0.325      0.745      -0.033       0.024\n",
      "sex           -0.0185      0.031     -0.590      0.555      -0.080       0.043\n",
      "age            0.2119      0.014     15.327      0.000       0.185       0.239\n",
      "urban          0.0703      0.016      4.393      0.000       0.039       0.102\n",
      "year_2004     -0.1884      0.085     -2.206      0.027      -0.356      -0.021\n",
      "year_2005     -0.2549      0.082     -3.125      0.002      -0.415      -0.095\n",
      "year_2006     -0.2353      0.081     -2.907      0.004      -0.394      -0.077\n",
      "year_2007     -0.4141      0.084     -4.921      0.000      -0.579      -0.249\n",
      "year_2008     -0.2392      0.082     -2.927      0.003      -0.399      -0.079\n",
      "year_2009     -0.3639      0.082     -4.446      0.000      -0.524      -0.204\n",
      "year_2010     -0.4284      0.085     -5.068      0.000      -0.594      -0.263\n",
      "year_2011     -0.3852      0.084     -4.564      0.000      -0.551      -0.220\n",
      "year_2012     -0.3844      0.084     -4.551      0.000      -0.550      -0.219\n",
      "year_2013     -0.3789      0.086     -4.413      0.000      -0.547      -0.211\n",
      "year_2014     -0.3828      0.084     -4.557      0.000      -0.547      -0.218\n",
      "year_2016     -0.5939      0.094     -6.298      0.000      -0.779      -0.409\n",
      "year_2018     -1.0039      0.102     -9.856      0.000      -1.204      -0.804\n",
      "year_2021      0.0004      0.082      0.004      0.996      -0.161       0.162\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "# ---------- 2) panel prep (two-way FE) ----------\n",
    "panel_numeric_cols = [\"id\",\"year\",\"ideology\",\"divorce\",\"housing_price\",\"income_1\",\n",
    "                      \"homeownership\",\"education\",\"sex\",\"age\",\"religion\"]\n",
    "for c in panel_numeric_cols:\n",
    "    if c in panel_df.columns:\n",
    "        panel_df[c] = pd.to_numeric(panel_df[c], errors=\"coerce\")\n",
    "\n",
    "# safe logs\n",
    "panel_df[\"ln_house_price\"] = safe_log(panel_df[\"housing_price\"])\n",
    "panel_df[\"ln_income\"]      = safe_log(panel_df[\"income_1\"])\n",
    "\n",
    "# key vars must be present\n",
    "panel_df = panel_df.dropna(subset=[\"id\",\"year\",\"ideology\",\"divorce\"])\n",
    "panel_df = panel_df.set_index([\"id\",\"year\"]).sort_index()\n",
    "\n",
    "# design matrix: DO NOT add constant in FE\n",
    "X1_cols = [\"divorce\",\"ln_income\",\"homeownership\",\"ln_house_price\",\n",
    "           \"education\",\"sex\",\"age\",\"religion\"]\n",
    "X1 = panel_df[X1_cols]\n",
    "y1 = panel_df[\"ideology\"]\n",
    "\n",
    "# fit two-way FE\n",
    "res1 = PanelOLS(y1, X1, entity_effects=True, time_effects=True,\n",
    "                drop_absorbed=True).fit(cov_type=\"robust\")\n",
    "\n",
    "# full results for FE model\n",
    "print(\"=== Two-Way Fixed Effects Model Results ===\")\n",
    "print(res1.summary)\n",
    "\n",
    "# ---------- 3) cross-sectional prep (Logit) ----------\n",
    "cross_numeric_cols = [\"conservative\",\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"]\n",
    "for c in cross_numeric_cols:\n",
    "    if c in cross_df.columns:\n",
    "        cross_df[c] = pd.to_numeric(cross_df[c], errors=\"coerce\")\n",
    "\n",
    "cross_df = cross_df.dropna(subset=[\"conservative\",\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"])\n",
    "\n",
    "y2 = cross_df[\"conservative\"].astype(int)\n",
    "\n",
    "# year dummies; ensure float dtype for statsmodels\n",
    "X2 = pd.get_dummies(cross_df[[\"divorce\",\"income\",\"education\",\"sex\",\"age\",\"urban\",\"year\"]],\n",
    "                    columns=[\"year\"], drop_first=True)\n",
    "X2 = X2.apply(pd.to_numeric, errors=\"coerce\").astype(float)\n",
    "X2 = sm.add_constant(X2, has_constant=\"add\")\n",
    "\n",
    "# Logit model\n",
    "logit_res = sm.Logit(y2, X2).fit()\n",
    "\n",
    "# full results for Logit model\n",
    "print(\"\\n=== Logistic Regression Model Results ===\")\n",
    "print(logit_res.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05a8acc5-c02d-4c0b-8f22-fb67419d1f96",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7b2158ac-aad6-431f-b71b-7f4b744e594c",
   "metadata": {},
   "source": [
    "# 3. Analysis 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "67a39968-16de-415e-9405-80b697cbc3cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"\\df_analysis_3.xlsx\")   # 복지패널"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "da9a9813-c578-441f-8535-8488d3b9899e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== District FE results (robust SEs) ===\n",
      "                           PanelOLS Estimation Summary                            \n",
      "==================================================================================\n",
      "Dep. Variable:     conservativeelectors   R-squared:                        0.5858\n",
      "Estimator:                     PanelOLS   R-squared (Between):             -62.721\n",
      "No. Observations:                  1260   R-squared (Within):               0.5858\n",
      "Date:                  Wed, Sep 17 2025   R-squared (Overall):             -59.515\n",
      "Time:                          16:14:57   Log-likelihood                    1989.6\n",
      "Cov. Estimator:                  Robust                                           \n",
      "                                          F-statistic:                      140.99\n",
      "Entities:                           253   P-value                           0.0000\n",
      "Avg Obs:                         4.9802   Distribution:                  F(10,997)\n",
      "Min Obs:                         1.0000                                           \n",
      "Max Obs:                         5.0000   F-statistic (robust):             154.59\n",
      "                                          P-value                           0.0000\n",
      "Time periods:                         5   Distribution:                  F(10,997)\n",
      "Avg Obs:                         252.00                                           \n",
      "Min Obs:                         252.00                                           \n",
      "Max Obs:                         252.00                                           \n",
      "                                                                                  \n",
      "                              Parameter Estimates                               \n",
      "================================================================================\n",
      "              Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI\n",
      "--------------------------------------------------------------------------------\n",
      "divorce_rate    -0.0464     0.0092    -5.0647     0.0000     -0.0643     -0.0284\n",
      "ln_economy       0.2521     0.0463     5.4437     0.0000      0.1612      0.3429\n",
      "palma           -0.0603     0.0113    -5.3270     0.0000     -0.0825     -0.0381\n",
      "high_edu         0.0044     0.0002     20.137     0.0000      0.0040      0.0049\n",
      "house           -0.0087     0.0042    -2.0609     0.0396     -0.0170     -0.0004\n",
      "mean_age        -0.0001     0.0024    -0.0501     0.9600     -0.0048      0.0046\n",
      "gender_ratio    -0.0004   9.05e-05    -4.1626     0.0000     -0.0006     -0.0002\n",
      "urbanization    -0.0006     0.0003    -2.0794     0.0378     -0.0012  -3.405e-05\n",
      "turn_out         0.4670     0.0927     5.0395     0.0000      0.2852      0.6489\n",
      "presidential     0.0034     0.0111     0.3041     0.7611     -0.0184      0.0251\n",
      "================================================================================\n",
      "\n",
      "F-test for Poolability: 23.890\n",
      "P-value: 0.0000\n",
      "Distribution: F(252,997)\n",
      "\n",
      "Included effects: Entity\n"
     ]
    }
   ],
   "source": [
    "DIV_GRID  = np.linspace(0.5, 3.5, 31)  # divorce_rate from 0.5 to 3.5\n",
    "PRES_SCENARIOS = [0]  # set to [0,1] if you want separate lines for non-presidential vs preside\n",
    "\n",
    "# Stata: gen ln_economy = ln(average_premium)\n",
    "df[\"ln_economy\"] = np.where(df[\"average_premium\"] > 0,\n",
    "                            np.log(df[\"average_premium\"]), np.nan)\n",
    "\n",
    "# Stata: encode district, gen(district_1)\n",
    "# We'll factorize district into an integer id\n",
    "df[\"district_1\"], _ = pd.factorize(df[\"district\"], sort=True)\n",
    "\n",
    "# Minimal cleaning for variables used\n",
    "needed = [\"district_1\",\"year\",\"conservativeelectors\",\"divorce_rate\",\"ln_economy\",\n",
    "          \"palma\",\"high_edu\",\"house\",\"mean_age\",\"gender_ratio\",\n",
    "          \"urbanization\",\"turn_out\",\"presidential\"]\n",
    "for c in needed:\n",
    "    df[c] = pd.to_numeric(df[c], errors=\"coerce\")\n",
    "df = df.dropna(subset=needed)\n",
    "\n",
    "# Panel index (district FE)\n",
    "df = df.set_index([\"district_1\",\"year\"]).sort_index()\n",
    "\n",
    "# =========================================================\n",
    "# 2) Build design matrix (FE, no constant)\n",
    "# Stata model:\n",
    "# xtreg conservativeelectors divorce_rate ln_economy palma high_edu house\n",
    "#       mean_age gender_ratio urbanization turn_out i.presidential, fe\n",
    "# =========================================================\n",
    "X_cols = [\"divorce_rate\", \"ln_economy\", \"palma\", \"high_edu\", \"house\",\n",
    "          \"mean_age\", \"gender_ratio\", \"urbanization\", \"turn_out\", \"presidential\"]\n",
    "\n",
    "y = df[\"conservativeelectors\"]\n",
    "X = df[X_cols]\n",
    "\n",
    "# Fit FE over districts (entity_effects=True), time effects not included (matches Stata spec)\n",
    "fe_res = PanelOLS(y, X, entity_effects=True, time_effects=False, drop_absorbed=True)\\\n",
    "           .fit(cov_type=\"robust\")\n",
    "\n",
    "print(\"\\n=== District FE results (robust SEs) ===\")\n",
    "print(fe_res.summary)\n",
    "\n",
    "# Extract params and robust covariance\n",
    "beta = fe_res.params\n",
    "V    = fe_res.cov\n",
    "\n",
    "# =========================================================\n",
    "# 3) Adjusted predictions as divorce_rate varies (0.5 → 3.5)\n",
    "#    Hold other regressors at their sample means.\n",
    "#    We produce either one line (PRES_SCENARIOS=[0]) or two lines (0 and 1).\n",
    "# =========================================================\n",
    "# Means of covariates (excluding divorce_rate and presidential which we set explicitly)\n",
    "means = X.drop(columns=[\"divorce_rate\",\"presidential\"]).mean()\n",
    "\n",
    "def predict_divorce_curve(div_vals, pres_val):\n",
    "    \"\"\"\n",
    "    Returns predicted conservative vote share and 95% CI as divorce_rate varies,\n",
    "    holding other covariates at their means and presidential fixed at pres_val.\n",
    "    \"\"\"\n",
    "    # Compose design rows for each divorce value\n",
    "    # order must match X_cols\n",
    "    base_vec = {**means.to_dict(), \"presidential\": pres_val}\n",
    "    # Ensure all keys present\n",
    "    for k in X_cols:\n",
    "        if k not in base_vec and k not in [\"divorce_rate\"]:\n",
    "            base_vec[k] = 0.0\n",
    "\n",
    "    preds, se_pred = [], []\n",
    "    # Build Jacobian for linear prediction: grad = x (since linear in parameters)\n",
    "    for dval in div_vals:\n",
    "        x_row = np.array([\n",
    "            dval,\n",
    "            base_vec[\"ln_economy\"],\n",
    "            base_vec[\"palma\"],\n",
    "            base_vec[\"high_edu\"],\n",
    "            base_vec[\"house\"],\n",
    "            base_vec[\"mean_age\"],\n",
    "            base_vec[\"gender_ratio\"],\n",
    "            base_vec[\"urbanization\"],\n",
    "            base_vec[\"turn_out\"],\n",
    "            base_vec[\"presidential\"]\n",
    "        ], dtype=float)\n",
    "\n",
    "        # Predicted \"within\" component: x_row' * beta (entity FE are absorbed; we report adjusted predictions)\n",
    "        yhat = float(np.dot(x_row, beta.values))\n",
    "        preds.append(yhat)\n",
    "\n",
    "        # Delta method for SE of linear prediction: SE = sqrt(x' V x)\n",
    "        se = float(np.sqrt(np.dot(x_row, np.dot(V.values, x_row))))\n",
    "        se_pred.append(se)\n",
    "\n",
    "    preds = np.array(preds)\n",
    "    se_pred = np.array(se_pred)\n",
    "    ci_l = preds - 1.96 * se_pred\n",
    "    ci_u = preds + 1.96 * se_pred\n",
    "    return preds, ci_l, ci_u\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c39d255-a4d0-4392-871b-86bc3a7149b2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb703928-f3ad-4730-b219-0bd2bb3240c0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
